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FOREWORD 

This  is  one  of  Che  two  Final  Reports  of  a  research  sponsored  by  the  Air 
Force  Office  of  Scientific  Research,  and  it  corresponds  to  the  analytical  part 
of  the  investigation.  The  second  Final  Report,  corresponding  to  the  experi¬ 
mental  part,  is  being  issued  separately  by  the  Geotechnical  Engineering  Center 
of  the  University  of  Texas  at  Austin,  as  follows: 

"Investigation  of  Low-Amplitude  Compression  Wave  Velocity  in 
Anisotropic  Material, ”  by  Shannon  H.H.  Lee  and  Kenneth  H.  Stokoe,  II, 
Report  GR86-6,  August  1986,  320  pages. 
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ABSTRACT 


/ 

!  -  / 

The  need  for  a  micromechanical  approach  to  modelling  the  stress-s train 
response  of  granular  soil  is  discussed  and  justified.  The  report  focuses  on 

t 

the  small  shear  strain  (y  _<  0.01%)  behavior,  and  investigates  the  validity  of 

V 

analytically  modelling  uniform,  rounded-grained  quartz  sands  by  arrays  of 

rr,>  . 

identical  elastic  quartz  spheres .  -zr— first  the  stress-strain  proper- 

ties  ofy^ix-  regular  arrays  of  spheres  are  studied^in--scnne  detatij  with  focus 
on  isotropic  and  transversely  isotropic  boundary  loading. 

An  analytical  procedure  is  established  for  determining  the  elastic  moduli 
of  a  random  assemblage  of  equal  elastic  spheres  of  arbitrary  mean  porosity, 
subjected  to  isotropic  confining  pressure.  The  procedure  uses  the  properties 
of  the  regular  arrays  already  described, 'd.t  accounts  for  the  spatial  distri¬ 
bution  of  porosity,  and  <it  calculates  the  macroscopic  moduli  through  the  ^Self 
/ 

Consistent  Method.  The  procedure  was  applied  to  compute  the  shear  and  bulk 
moduli  of  assemblages  of  quartz  spheres*  which  were  then  compared  with  static 
and  dynamic  measurements  on  quartz  sands  from  the  literature.  The  theoretical 
sands  are  significantly  stiffer  than  the  actual  soils  due  to  the  lower  number 
of  effective  contacts  in  actual  sands.  However,  excellent  agreement  was  found 
with  resonant  column  shear  modulus  measurements  on  Ottawa  sand,  after  subject¬ 
ing  it  to  a  large  number  of  cycles  of  shear  prestraining,  which  increased  the 
number  of  contacts  toward  the  theoretical  value. 
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Section  1 
INTRODUCTION 


The  main  objective  of  this  report  is  to  present  a  simple,  yet  rigorous, 
particulate  mechanics  model  of  the  stress-strain  response  of  granular  soil 
under  isotropic  boundary  loading  at  very  small  shear  strains,  y,  of  the  order 
of  10-^%  or  less.  The  proposed  model  idealizes  sand  as  a  combination  of  regu¬ 
lar  arrays  of  elastic,  rough  spheres  and  uses  Mindlin's  formulation  for  the 
contacts.  In  turn,  this  is  the  first  stage  of  a  long-term  attempt  to  model 
sands  as  3-D  spatial  arrangements  of  regular  arrays  at  both  small  strains  (y 
Yt  =  0.01%)  and  large  strains  (y  >  yt  =  0.01%),  where  =  0.01%  is  the  thres¬ 
hold  strain  for  densif ication  and  pore  water  pressure  buildup^*).  A  second 
stage  will  include  anisotropically  loaded  granular  media,  and  the  ultimate  goal 
is  to  perform  2-D  and  3-D  computer  simulations  of  arrays  of  spheres  at  differ¬ 
ent  small  and  large  strain  ranges,  including  analytical  modelling  of  densif ica¬ 
tion  under  boundary  cyclic  loading. 

This  is  a  final  report  of  a  research  on  the  subject  performed  by  the 
authors  in  cooperation  with  a  U.  of  Texas  team  headed  by  Prof.  Stokoe. 

"Elastic  constants"  of  interest  at  very  small  strains  include  the  shear 
and  bulk  moduli  and  the  Poisson's  Ratio(s).  Experimental  results  and  basic 
considerations  indicate  that  these  "constants"  depend  on  both  the  void  ratio  of 
the  soil  and  the  state  of  confining  stresses.  The  variations  of  these  moduli 
and  of  the  damping  of  the  soil  with  an  applied  shear  strain  up  to  the  threshold 
are  also  of  interest,  as  is  the  value  of  the  threshold  strain  itself  at  which 

The  concepts  of  small  and  large  shear  strains  as  used  here  are  consistent 
with  usual  Soil  Dynamics  terminology,  but  they  do  not  coincide  with  that  of 
traditional  Soil  Mechanics,  where  much  greater  strains,  about  1%  or  larger 
are  usually  of  interest. 
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gross  sliding  occurs  at  the  grains'  contacts. 

These  small  strain  soil  parameters  are  very  important  in  geotechnical 
engineering  problems  involving  cyclic  loading  or  wave  propagation  in  the  soil, 
such  as:  ocean  wave  loading,  soil  structure  interaction,  site  response, 
ground  settlement  and  liquefaction  during  earthquakes.  Due  to  this,  a  great 
number  of  experimental  studies  of  small  strain  behavior  have  been  performed, 
and  correlations  have  been  developed  for  practical  use.  Especially  important 
are  the  equations  for  the  shear  modulus  at  very  small  strains,  Gmax,  in  sands 
developed  by  Hardin  and  Richart  (1963)  and  Seed  and  Idriss  (1970)  on  the 
assumption  that  these  soils  can  be  treated  as  elastic  isotropic  solids.  In 
both  correlations,  discussed  in  more  detail  in  Section  2,  Gmax  =  A*(o0)®*^, 
with  oj,  02,  03  being  the  effective  principal  stresses,  o0  =  (oj  +  02  +  a3)/3 
is  the  mean  effective  stress  and  A  is  a  soil  constant  which  depends  on  void 
ratio  or  relative  density.  Both  correlations  assume  that  Gmax  (and  thus,  also 
the  shear  wave  velocity,  Vg  =  (Gmax/p)^2),  is  the  same  for  isotropically  or 
anisotropically  loaded  sand,  provided  that  the  mean  stress  oQ  is  the  same; 
also,  both  correlations  assume  that  for  the  anisotropic  loading  case  Gmax  and 
Vg  do  not  change  with  direction. 

These  assumptions  for  Gmax  in  sands  have  been  challenged  more  recently  by 
the  experimental  results  obtained  by  Roesler  (1979),  Knox  et  al.  (1982)  and  Yu 
and  Richart  (1984),  as  discussed  in  Section  2.  Therefore,  a  main  motivation 
for  this  work  was  the  need,  suggested  by  those  experimental  findings,  for  a 
fresh  approach  to  our  basic  understanding  of  Gmax  and  other  small-strain  soil 
parameters.  Some  preliminary  analytical  results  previously  obtained  by  the 
senior  author  (Dobry  et  al. ,  1982)  had  shown  that  a  particulate  mechanics 
approach  was  very  well  suited  to  this  purpose,  and  should  be  the  basis  of  this 
fresh  approach. 


The  large  strain  (0.012  <  y  <  1%)  behavior  of  granular  soils  is  also  very 
important  in  engineering  problems  involving  cyclic  loading,  and  especially 
those  related  to  earthquakes.  At  these  strains,  the  stress-strain  behavior 
becomes  strongly  nonlinear  and  hysteretic,  and  rearrangement  of  particles  take 
place  producing  phenomena  such  as  densif ication  and  pore  water  pressure  build¬ 
up  (Silver  and  Seed,  1971;  Youd,  1972;  Dobry  et  al. ,  1982;  National  Research 
Council,  1985).  A  number  of  continuum  mechanics  models,  based  mostly  on  the 
Incremental  Theory  of  Plasticity,  try  to  simulate  this  behavior  and  have  been 
proposed  for  soils,  as  discussed  in  Section  3. 

A  summary  literature  review  of  previous  relevant  particulate  mechanics 
studies  is  presented  in  Section  4.  Many  of  these  past  investigations  have 
focused  on  the  very  large  strain  (y  >  1%)  and  strength  behavior  of  granular 
soils;  at  those  very  large  strains,  gross  sliding  and  rolling  of  the  grains  are 
main  contributors  to  the  overall  strain,  while  the  elasticity  of  the  particles 
and  contacts  play  a  minor  or  negligible  role.  On  the  other  hand,  for  the  small 
to  large  strain  ranges  of  interest  of  the  proposed  research,  the  elasticity  of 
the  particles  and  the  details  of  the  force-displacement  response  at  the  contacts 
are  very  significant  factors.  The  discussion  in  Section  4  includes  a  general 
model  recently  proposed  for  the  force-displacement  response  at  the  contact 
between  two  identical  elastic  spheres  (Seridi  and  Dobry,  1984). 

The  results  of  the  present  research  are  discussed  in  Sections  5  and  6. 
Section  5  presents  a  detailed  study  of  the  differential  stress-strain  relations 
for  various  regular  arrays  of  spheres.  Section  6  describes  an  application  of 
these  findings,  using  the  Self-Consistent  Method,  to  a  random  arrangement  of 
regular  arrays  subjected  to  isotropic  boundary  loading,  and  with  the  arrangement 
having  an  arbitrary  macroscopic  void  ratio. 
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LABORATORY  MEASUREMENTS  ON  SANDS  AT  SMALL  STRAINS 

Starting  around  1960,  a  number  of  cyclic  and  dynamic  laboratory  measure¬ 
ments  have  been  performed  to  determine  the  stress-strain  behavior  of  granular 
arrays  and  of  natural  sands  at  small  strains.  Properties  studied  have  in¬ 
cluded:  i)  maximum  shear  modulus  at  very  small  strains,  Gmax;  ii)  the  varia¬ 

tion  of  secant  modulus,  G,  with  shear  strain,  y;  iii)  the  Poisson's  ratio  of 
the  soil;  iv)  the  variation  of  shear  damping  ratio  with  strain;  and  v)  the 
threshold  shear  strain,  yt,  at  which  densif ication  and  pore  pressure  buildup 
start.  Many  of  these  tests  have  been  conducted  in  a  triaxial  cell,  on  sand 
specimens  isotropically  or  anisotropically  consolidated  under  a  biaxial  stress 
state  (02  =  03  or  02  =  03),  with  the  small  strain  measurements  performed  using 
the  pulse  method,  the  resonant  column  or  cyclic  torsional  techniques,  and  with 
particular  emphasis  on  shear  modulus  determinations.  Important  results  and 
state-of-the-art  summaries  of  these  modulus  measurements  have  been  presented 
by  Duffy  and  Mindlin  (1957),  Hardin  and  Richart  (1963),  Lawrence  (1965), 
Richart  et  al.  (1970),  Seed  and  Idriss  (1970),  Hardin  and  Drnevich  (1972), 
Woods  (1978),  Iwasaki  et  al.  (1978),  and  Tatsuoka  et  al.  (1979). 

As  mentioned  before,  "small  strains"  are  defined  here  by  the  condition  y 
<  yt,  as  in  this  range  the  original  geometry  of  the  granular  array  or  sand 
remains  essentially  unchanged,  with  very  few  or  no  particles  experiencing 
gross  sliding  or  rolling,  and,  thus,  with  the  macroscopic  strain  of  the  array 
being  controlled  by  the  elastic  deformations  of  the  particles  and  by  localized 
slips  within  the  areas  of  contact  areas  between  particles.  In  many  sands,  yt 
=  10-^  =  10~2%;  measurements  and  studies  of  y^  have  been  presented  by  Drnevich 

and  Richart  (1970),  Youd  (1972),  Pyke  (1973),  Dobry  et  al.  (1980,  1981,  1981a, 
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1982),  Dyvik  et  al.  (1982),  Oner  (1984)  and  National  Research  Council  (1985). 

On  the  basis  of  laboratory  measurements,  Hardin  and  Richart  (1963) 
proposed  the  following  expression  for  Gmax: 

Gmax  =  f(e)(a0)°*5  (1) 

where  e  =  void  ratio,  f(e)  =  2630(2. 1 7-e)^/(l+e)  for  round-grained  sands,  and 
f(e)  =  1230(2. 97-e)2/(i+e)  for  angular-grained  sands;  both  Gmax  and  oG  are  in 
psi  in  Eq.  1. 

Seed  and  Idriss  (1970)  proposed  the  alternative  expression: 

Gmax  =  K2max  (o0)G*^  (2) 

where  Gmax  and  a0  are  in  psf,  and  K2max  a  function  of  relative  density. 

Eqs.  1-2  reflect  the  conclusion  of  these  and  other  studies,  that  Gmax  and 
Vg  are  mainly  a  function  of  void  ratio  or  relative  density,  and  of  the  mean 
effective  normal  stress  oQ.  Other  variables,  such  as  static  shear  stresses, 
stress  history,  stress  path  (compression  versus  extension  loading),  frequency 
of  cyclic  loading,  degree  of  saturation,  were  found  to  have,  either  a  small 
effect  or  no  effect  at  all  (Richart  et  al.,  1970,  Yu  and  Richart,  1984).  One 
exception  is  that  a  large  number  of  shear  prestraining  cycles  at  strains  larger 
than  the  threshold  was  found  to  increase  Gmax  significantly  (Drnevich  and 
Richart,  1970). 

It  is  useful  to  make  explicit  some  of  the  implications  of  Eqs.  1-2  for  an- 
isotropically  loaded  dry  sand,  either  for  the  general  "triaxial"  case  in  which 
ai  *  02  *  03  or,  as  is  very  usually  the  case  in  the  field,  or  for  the  ’'biaxial' 
case,  oj  *0  2  *  03.  (The  bars  have  now  been  dropped  from  the  stresses,  as  for 
a  dry  sand,  o  =  a).  These  implications  are: 


i)  Gmax  and  Vs  depend  equally  on  ai,  02  and  °3* 
li)  For  a  shear  wave  propagating  in  the  sand,  the  value  of  Vs  is  the  same 
whatever  the  directions  of  propagation  and  polarization  of  the  wave. 

Implication  ii)  is  equivalent  to  assume  that,  at  very  small  strains,  the 
anisotropically  loaded  sand  can  be  modeled  as  an  isotropically  elastic  mater¬ 
ial,  defined  by  the  two  elastic  constants  Gmax  and  Poisson's  Ratio,  v.  Under 
the  usual  additional  assumption  that  for  a  dry  sand,  v  =  0.3  to  0.4  is  a 
constant  independent  of  confining  stresses,  a  set  of  implications  similar  to  i) 
and  ii)  can  be  obtained,  but  now  for  a  dilational,  P-wave,  travelling  in  this 
dry  sand.  If  we  call  D,  Vp  the  constrained  modulus  and  P-wave  velocity,  the 
corresponding  implications  are: 

iii)  D  and  Vp  depend  equally  on  oj,  02  an^  °3»  with  0  a(aQ)^*^  and  Vp 
a(o0)0*25. 

iv)  The  value  of  Vp  is  the  same  whatever  the  direction  of  propagation  of 
the  P-wave. 

Five  recent  experimental  studies  have  attempted  to  verify  in  detail  this 
formulation  by  Hardin/Ri chart  and  Seed/Idriss,  and  specifically  the  validity 
of  implications  i)  through  iv)  above  for  anisotropically  loaded  dry  sand. 

Schmertmann  (1978)  measured  Vp  and  Vg  in  several  directions  in  a  large  dry 
sand  specimen  (4  ft  diameter  by  4  ft  high).  In  these  tests,  a  biaxial  state  of 
stress  could  be  achieved,  =  ov  t  02  =  03  =  o^,  with  av,  o^  =  vertical,  hori¬ 
zontal  stresses,  with  the  stresses  varying  between  5  to  20  psi,  and  with  a 
stress  ratio,  0^/03  =  1  to  3.  He  found  that  there  was  a  slight  amount  of  inher¬ 
ent  anisotropy  (different  wave  velocities  in  the  horizontal  and  vertical  direc¬ 
tions  when  av  =  a^).  He  also  found  that  for  constant  oQ  =  l/3(ov  +  2  o^)  and 
variable  0^/03,  Vg  varied  less  than  10%,  thus  verifying  the  basic  Hardin/Richart 
assumption  as  a  first  approximation  for  Vg  in  this  biaxial  case.  However,  Vp 


was  strongly  affected  by  01/03.  The  results  suggested  that,  for  P-waves  propa¬ 
gating  in  the  vertical  direction,  Vp  depended  more  on  ov  than  on  o0. 

Roesler  (1979)  measured  Vg  using  a  1  ft^  cubical  dry  sand  sample.  In 
these  tests,  a  true  triaxial  state  of  stresses  was  achieved.  Test  pressures 
ranged  from  5.8  psi  to  23  psi,  with  01/03  =  1  to  1.8.  He  propagated  the 
shear  waves  along  either  of  the  principal  stress  directions  (oa),  with  particle 
motions  polarized  in  another  principal  direction  (o^).  The  third  principal 
direction,  or  out-of-plane  direction,  is  neither  a  direction  of  propagation  nor 
polarization  (oc).  Roesler  found  that  his  results  followed  the  law: 


B  o  0-1*9  0.  0.107  0 


where  B  =  constant.  These  results  are  illustrated  in  Fig.  1.  For  isotropic 
confinement  (a0=aa=ai,=  ac)  they  do  confirm  the  Hardin-Ri chart  law  that  Gmax 
<*(a0)0*5  and  Vs  a (a0)®*^,  as  1.49  +  0.107  =  0.256.  However,  for  the  general 
case,  Eq.  3  contradicts  Eqs.  1-2,  in  that  now  Vg  is  completely  independent  of 
ac.  Also,  Roesler's  results  for  this  case  indicate  that  Vs  is  a  function  of 
direction  and  the  sand  cannot  be  treated  as  an  isotropic  elastic  material; 
therefore,  more  than  two  elastic  constants  are  necessary  to  define  it. 

Stokoe  et  al.  (1980)  developed  at  the  University  of  Texas  at  Austin 
(U.T.)  a  large  scale,  7  ft-*  cubical  triaxial  facility,  for  the  specific  pur¬ 
pose  of  measuring  Vg  and  Vp  in  dry  sand.  In  this  facility,  a  triaxial  state 
of  stress,  a\  *  03  *  03  can  also  be  achieved.  All  tests  performed  to  date 
at  U.T.  have  used  a  local  medium  to  fine,  washed  mortar  sand  classified  as  SP, 
with  effective  grain  size,  Djq  =  0.28  mm  and  a  uniformity  coefficient  Cu  =  1.7. 


The  sand  is  placed  by  the  raining  technique  and  is  tested  dry.  Values  of 
principal  stresses  used  have  ranged  between  10  psi  and  40  psi,  with  the  stress 
ratio,  0^/0 3  =  1  to  4.  Knox  et  al.  (1982)  used  this  facility  to  study  Vg  and. 


similarly  to  Roesler,  they  propagated  the  shear  waves  along  one  of  the 
principal  stresses  (aa),  and  polarized  parallel  to  another  principal  stress 
(oj,)*  They  expressed  Vs  as: 

Vs  =  F  oama  obmb  ocrac  (4) 


where  F  =  constant.  They  found  values  of  ma  =  mb  =  0.09  to  0.12,  and  me  =  0  to 
0.01.  Except  for  some  minor  differences,  these  results  are  identical  to  those 
of  Roesler,  including  the  independence  of  Vg  on  ac.  Kopperman  et  al.  (1982) 
used  the  same  U.T.  facility  and  sand  to  study  P-waves  and  concluded  that: 


VP  = 


L  a 


0.22 


(5) 


where  L  =  constant.  The  insensitivity  of  Vp  to  variations  in  the  stresses 
and  oc  perpendicular  to  wave  propagation  is  illustrated  by  Fig.  2.  Again,  and 
similar  to  Schmertmann's  findings,  these  results  indicate  that  Vp  is  strongly 
dependent  on  direction  of  propagation  when  the  sand  is  anisotropically  loaded. 

Yu  and  Richart  (1984)  performed  resonant  column  tests  on  three  sands  sub¬ 
jected  to  a  biaxial  state  of  stress.  Their  results  essentially  agreed  with 
those  of  Roesler  and  Knox;  however,  they  found  some  effect  of  the  stress  ratio 
on  the  results.  They  proposed  for  Gmax  the  expression: 


G  =  CP  0.49 
umax  ora 


av° * 2 6  0h0'25  U“a  Kn2> 


(6) 


where  Oconstant,  Pa=atmospheric  pressure,  a=0.15  to  0.23,  with  a  mean  value  of 
0*18.  ^=(01/03- * ) / [ (al/a3)max-l  1 »  anc*  ^al/°2^max  corresponds  to  shear  failure 
of  the  sand.  Except  for  the  factor  1-a  Kn2,  which  is  usually  between  0.8  and 
1,  Eq .  6  is  consistent  with  Eqs.  3  and  4  proposed  by  Roesler  and  Knox. 


Therefore,  all  of  these  results  clearly  indicate  that  implications  i) 
through  iv)  above,  associated  with  the  currently  used  correlations  for  V  and 


Vs  in  sands,  need  Co  be  revised  and  upgraded.  In  most  cases  of  practical 
interest,  sands  are  anisotropically  loaded,  and  thus  more  than  two  elastic 
constants  may  be  needed  to  specify  the  behavior  of  the  soil  at  very  small 
strains.  For  the  typical  biaxial  state  of  stresses  existing  in  the  field  under 
geostatic  conditions,  for  which  all  horizontal  normal  stresses  are  equal  but 
different  from  the  vertical  stress,  the  sand  will  behave  as  a  cross-anisotropic 
elastic  solid  and  5  elastic  constants  will  be  generally  needed  (Love,  1927, 
Sokolnikoff,  1956).  In  the  more  general  case  of  *  02  *  a 3,  as  it  may  happen 
in  the  soil  under  a  structure,  the  sand  can  be  described  as  an  orthotropic 
elastic  solid,  with  three  planes  of  elastic  symmetry,  and  a  total  of  9  elastic 
constants  are  needed  (Sokolnikoff,  1956). 

The  previous  discussion  focused  on  the  elastic  properties  of  sand  at  very 
small  strains,  y  -  10“^%,  and  especially  on  Vp  and  Vs  measurements.  If  larger 
loads  and  strains  are  applied  to  a  dry  granular  soil,  compression-wave  type 
loading  induces  a  nonlinear  locking  stress-strain  response,  while  shear-wave 
type  loading  induces  a  yielding  response  (see  Fig.  3).  This  behavior  is 
obviously  associated  with  the  particulate  nature  of  the  soil  (Seed  and  Idriss, 
1970;  Hardin  and  Drenvich,  1972).  During  cyclic  shear  loading  in  sand,  stress- 
strain  hysteresis  loops  are  generated  such  as  shown  in  Fig.  4;  these  loops  are 
essentially  strain-rate  and  frequency  independent.  For  small  strains,  y  <  yt  = 
10  the  hysteretic  loop  repeats  itself  cycle  after  cycle,  and  no  permanent 
volumetric  strain  is  observed,  thus  suggesting  an  essentially  non-destructive 
though  nonlinear  behavior,  controlled  mainly  by  the  response  of  the  contacts 
between  the  grains,  and  with  no  coupling  between  shear  and  volumetric  strains. 
At  shear  strains,  y  >  Yt»  although  the  overall  behavior  remains  approximately 
the  same,  densif ication  occurs,  and  there  is  also  some  increase  in  stiffness, 
with  the  shear  stress-strain  curve  and  the  tips  of  the  loop  going  up  a  little 


as  the  number  of  cycles  increases  (see  Fig.  4).  For  the  strain  range 
0.01%  <  y  <0.1%  to  1%,  the  monotonic  and  cyclic  behavior  of  the  sand  is  always 
contractive,  that  is,  shear  strains  generate  exclusively  compressive  volumetric 
strains,  independently  of  the  density  of  the  sand.  However,  for  strains  larger 
than  about  y  =  1%,  a  mixture  of  contractive  and  dilative  behavior  is  measured 
in  dense  sands,  with  expansion  of  the  soil  occurring  during  part  of  the  cycle 
in  cyclic  shear  loading  (Youd,  1972).  At  all  strains,  and  for  both  shear  and 
compression-type  cyclic  loading,  the  stress-strain  response  of  dry  granular 
soil  is  strongly  dependent  on  the  level  of  normal  stresses  acting  on  it. 


Section  3 


STRESS-STRAIN  MATHEMATICAL  MODELLING 


A  significant  amount  of  research  has  been  directed  to  obtain  stress-strain 

constitutive  relations  for  cyclic  and  dynamic  loading  of  soil.  Most  of  these 

studies  have  modelled  the  soil  as  an  elastic-plastic  material,  using  as  a  basis 

tool  the  Incremental  Theory  of  Plasticity.  In  this  type  of  model,  which  is 

particularly  appropriate  for  dry  granular  soil,  the  total  strain  increment  is 

e  p 

equal  to  the  sum  of  elastic  and  plastic  strain  increments,  de  =  de  +  de  , 

ij  ij  ij 

with  all  de  being  strain  rate  independent  (Drucker  and  Prager,  1952;  Reyes, 

1966;  Chen,  1975;  Lade  and  Duncan,  1975;  Prevost,  1978;  Hardin,  1978).  Based 

on  the  Vp  measurements  by  Roesler  previously  described  in  Section  2,  Hardin 

e 

(1980)  suggested  the  following  expressions  for  de  in  dry  granular  soil: 

ij 


,  e  F(e)  *  ,  d<Jx  day  doz 

dex  - -  ( - v  - 1 - v  - J 
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e  e  e 

where  t  =  normal  strain,  r  *  2c  =  engineering  shear  strain,  and  four 
x  xv  xy 

additional  equations  are  obtained  by  permutation  of  subscripts.  In  these 
equations,  -j  x ,  jv  arid  j  7  =  normal  stresses;  t  ^  =  a  shear  stress;  Pa  = 
atmospheric  pressure;  and  F'  e  1  =  ‘.i  +  0./  e4- ,  where  e  =  void  ratio.  Eqs.  7 

contain  five  el  as*  ;  -os*  *uts  ’  >x,  Sv,  Sz ,  SXy  and  v );  based  on  Roesler's 
experimental  ,  >>wer  >T  stress  n  »  0.5  was  proposed. 

A  vari  et  v  tss  i  •  »•  :  arid  nuiiassoc  iated  flow  roles  have  been  proposed 

for  the  plastic  :  r  a :  ;  or  renn:  ,  of  the  form: 


where  X  =  coefficient  of  proportionality,  and  gCo^j)  is  the  plastic  potential 
function,  which  may  or  may  not  coincide  with  the  yield  function  fCajj)  at  which 
plastic  strains  develop.  Figure  5  shows  the  shapes  of  a  number  of  plastic 
potential  surfaces  proposed  for  soil  by  different  authors. 

In  the  simplest  type  of  elastic-plastic  model,  there  is  only  one  yield 
(failure)  surface.  For  stresses  below  that  surface,  the  behavior  is  assumed  to 
be  perfectly  elastic.  However,  granular  and  other  soils  develop  plastic  strains 
even  at  the  small  shear  strains  of  interest  to  this  report.  To  allow  for 
this  behavior,  a  wide  variety  of  strain-hardening  laws  have  been  proposed, 
including  families  of  yield  surfaces  and  specific  strain-hardening  yield  rules. 
In  some  of  these  models,  the  elastic  region  is  completely  eliminated,  thus 
allowing  for  plastic  flow  at  very  low  levels  of  stress  and  strain  (Mroz,  1967; 
Prevost,  1977).  One  of  the  earliest  developments  included  various  cap  models, 
based  on  the  work  done  at  Cambridge  University  by  Roscoe  and  his  co-workers 
(i.e.,  Roscoe,  1970).  This  includes  the  models  proposed  in  several  papers  by 
DiMaggio  and  Sandler  (1971),  which  have  been  widely  used  for  dynamic  analyses 
of  soil  response  to  explosions.  Several  capped  yield  models  are  included  in 
Fig.  5. 

An  important  aspect  of  the  development  of  elastic-plastic  models  Is  the 
definition  of  the  strain-hardening  law,  which  defines  the  modifications  of  the 
yield  surface(s)  in  the  course  of  the  plastic  flow.  This  is  especially 
critical  for  cyclic  loading,  where  the  type  of  strain-hardening  determines  the 
stress-strain  behavior  after  load  reversals.  In  most  of  the  models  described 
above,  which  were  originally  developed  for  monotonic  loading,  isotropic  strain¬ 
hardening  is  assumed  (Hill,  1950),  with  the  yield  surfaces  expanding  as  the 


stresses  increase  (Fig.  6).  When  isotropic  hardening  is  assumed,  a  large 
amount  of  load  reversal  is  required  for  additional  yielding  to  occur,  in  con¬ 
tradiction  with  the  observed  behavior  of  experimental  hysteresis  loops  such  as 
shown  in  Fig.  4. 

A  better  alternative  for  earthquake  loading  is  provided  by  the  kinematic 
strain-hardening  law,  sketched  in  Figure  7.  The  kinematic  model  was  originally 
proposed  by  Ishlinsky  (1954)  and  Prager  (1955).  Iwan  (1967)  proposed  a  reo- 
logical  representation  for  the  stress-strain  model,  constituted  by  infinite 
e lasto-plastic  elements,  placed  in  series  or  in  parallel.  This  model  is  a  non- 
frictional  one,  with  all  the  nested  yield  surfaces  being  circular  cylinders  in 
principal  stress  space.  Mroz  (1967,  1969)  proposed  a  general  model  for  elastic- 
plastic  materials,  composed  also  of  a  field  of  yield  surfaces,  with  a  combina¬ 
tion  of  kinematic  and  isotropic  strain-hardening  laws.  Prevost  (Prevost  and 
Hoeg,  1975;  Prevost,  1977,  1978),  Mroz,  Norris  and  Zienkiewicz  (1978,  1979)  and 
Vicente  and  Dobry  (1985),  have  proposed  to  use  this  model  to  predict  static  and 
cyclic  behavior  of  soils.  The  model  is  flexible  enough  to  allow  its  adaptation 
to  the  cases  of  drained  and  undrained  loading,  and  to  incorporate  important 
large  strain  cyclic  phenomena  such  as  densif ication,  liquefaction  and  stiffness 
degradation.  Anisotropically  loaded  soils  are  represented  by  nonsymraetric 
nested  surfaces  in  stress  space.  Under  cyclic  shear  loading,  the  strain-harden¬ 
ing  behavior  is  basically  kinematic  for  the  reasons  described  above.  A  simul¬ 
taneous  isotropic  hardening  (or  softening)  is  allowed  with  the  corresponding 
expansion  or  contraction  of  the  yield  surfaces  as  cyclic  loading  develops.  Thi: 
isotropic  expansion  (contraction)  thus  could  simulate  in  dry  granular  soil  the 
observed  increase  (f  stiffness  caused  by  cyclic  loading  above  yt  =  10~2%. 


Section  4 


THE  MICRCMECHANICAL  APPROACH 

Eqs.  1-2  for  Gmax  and  similar  expressions  for  other  small  strain  moduli  in 
dry  sands  assume  that  the  controlling  normal  stress  parameter  is  the  mean  stress 
o 0:  Gmax  =  This  functional  relationship,  selected  on  the  basis  of 

limited  experimental  evidence,  was  no  doubt  influenced  by  a  continuum  mechanics 
view  of  the  situation,  as  0o  is  proportional  to  the  first  invariant  of  the 
stress  tensor.  As  discussed  in  Section  2,  more  detailed  measurements  have 
revealed  the  elastic  anisotropy  of  a  dry  sand  subjected  to  different  principal 
stresses,  and  they  have  also  shown  that  the  functional  relationship  between 
principal  stresses  cfj,  02.  03.  on  one  hand,  and  Gmax  and  other  elastic  constants 
on  the  other,  is  not  Gmax  =  f(o0)  but  rather  Gmax  =  f(oa,aj,);  similarly,  for 
the  constrained  modulus,  D  =  f(aa).  Derivations  shown  later  in  Section  5.1 
for  a  simple  cubic  regular  array  of  elastic  rough  spheres  match  very  well  with 
those  recent  experimental  findings  in  sands.  This  strongly  suggests  that  a 
micromechanical  (particulate  mechanics)  approach  should  be  used  to  analytically 
simulate  and  generalize  the  experimental  observations. 

A  great  number  of  studies  have  been  performed  using  particulate  models  to 
understand  and  model  the  behavior  of  cohesionless  soils  and  other  granular 
materials.  Most  of  these  investigations  have  been  analytical,  but  they  have 
also  included  measurements  in  actual  granular  soils,  as  well  as  in  regular  or 
random  arrays  of  spheres  (3-D)  or  disks/rods  (2-D);  a  number  of  them  have 
dealt  with  the  load-deformation  characteristics  at  the  contact  between  two 
elastic  bodies  possessing  friction  (Mindlin's  problem).  State-of-the-art 
summaries  have  been  presented  by  Deresiewicz  (1958,  1973),  Mogami  (1965), 

Scott  and  Ko  (1969),  Richart  et  al.  (1970),  White  (1965),  Harr  (1977),  and 


Dobry  and  Grivas  (1978),  between  others.  The  proceedings  of  two  US/Japan 
seminars  on  the  mechanics  of  granular  materials  contain  excellent  papers  on 
the  subject  (Cowin  and  Satake,  1978;  Jenkins  and  Satake,  1983). 

A  number  of  these  studies  have  focused  on  the  probabilistic  aspects  and 
statistical  distributions  of  different  parameters  within  the  soil  or  granular 
medium,  and  their  effect  on  the  mechanical  behavior  of  the  array.  These  have 
included  investigations  of  the  orientations  of  the  individual  particles,  of 
the  spatial  distribution  of  porosity,  and  of  the  distribution  of  number, 
orientation  and  levels  of  force  transmitted  by  the  contacts,  conducted  between 
others  by  Smith  et  al.  (1929),  Dantu  (1957),  Field  (1963),  Mogami  (1965), 

Grivas  and  Harr  (1974),  Oda  (1974),  Yanagisawa  (1978),  Shahinpoor  and  Shahrpass 
(1982),  Nemat-Nasser  and  Mehrabadi  (1983),  and  Dobry  and  Petrakis  (1984). 

Many  of  those  analytical  investigations,  computer  simulations  and  obser¬ 
vations  have  focused  on  the  stress-strain  behavior  at  very  large  strains  and 
on  the  failure  of  dry  granular  media.  Because  of  this  very  large  strain 
nature  of  the  phenomena,  the  load-deformation  characteristics  of  the  particles 
contacts  have  played  a  minor  or  negligible  role,  and  the  emphasis  has  been  on 
changes  in  the  geometric  arrangement  of  the  grains  due  to  their  sliding  and 
rolling.  In  some  of  these  investigations  the  compliance  of  the  contacts  has 
been  eliminated  altogether  by  assuming  perfectly  rigid  particles.  Some  impor¬ 
tant  references  here  are  Rowe  (1962),  Morgenstern  (1963),  Horne  (1965), 

Konishi  (1978),  and  Oda  et  al.  (1983).  Cundall  and  Strack  (1983)  performed 
numerical  experiments  of  2-D  random  arrays  of  disks  using  an  explicit  finite 
difference  procedure  (see  Fig.  8).  In  these,  the  authors  successfully  simu¬ 
lated  "compression  triaxial”  loading  to  failure,  and  studied  in  detail  the 
spatial  distribution  of  contact  forces  and  the  distribution  and  relative  con¬ 
tributions  of  sliding  and  rolling  to  macroscopic  strain,  during  anisotropic 


deviatoric  loading  with  constant  03.  One  of  their  conclusions  for  this 
deviatoric  loading  is  that  the  major  principal  stress  oi  is  transmitted  mainly 
by  a  few  "stiff  chains"  of  particles  having  large  contact  forces,  with  the 
particles  in  between  chains  being  lightly  loaded;  sliding  and  rolling  occurs 
mainly  in  those  lightly  loaded  regions.  Oner  (1984)  worked  with  a  similar 
numerical  scheme  to  predict  the  observed  threshold  strain  Yc  at  which  sliding 
and  rolling  starts,  while  Dobry  et  al.  (1982)  used  a  regular  simple  cubic  array 
of  spheres  for  the  same  purpose. 

Of  special  interest  are  the  investigations  which  have  studied  the  stress- 

strain  behavior  of  granular  arrays  considering  the  elasticity  of  the  particles 

and  the  corresponding  compliances  at  the  contacts.  Most  of  these  studies  have 

assumed  spherical  grain  shapes  and  elastically  isotropic  grains  characterized 

by  three  material  constants  (two  elastic  and  one  friction  coefficient).  All 

of  these  investigations  have  used  the  normal  and  tangential  compliances  at  the 

contact  between  two  elastic  bodies,  derived  by  Hertz  (1882),  Cattaneo  (1938), 

and  Mindlin  and  his  co-workers  (Mindlin,  1949,  Mindlin  et  al.  1951,  Mindlin 

and  Deresiewicz,  1953).  Figures  9  and  10  show,  respectively,  the  distorsion  of 

two  spheres  subjected  to  normal  (N)  and  tangential  (T)  contact  loads,  and  the 

tangential  load-displacement  curve  for  constant  N.  As  noted  in  the  summary 

reviews  of  the  contact  theory  by  Deresiewicz  (1958,  1973)  and  Dobry  and  Grivas 

(1978),  Mindlin  and  his  co-workers  developed  the  basic  theoretical  framework 

of  the  contact  problem,  and  solved  it  for  some  special  force  time  histories; 

however,  the  general  problem  of  computing  the  displacements  for  a  contact 
+  -*■  + 

force  P  =  N  +  T,  whose  magnitude  and  direction  change  arbitrarily,  remained 
unsollved.  Only  very  recently,  Seridi  and  Dobry  (1984)  provided  a  general  and 
practical  solution  to  this  general  problem,  thus  making  it  possible  the  use  of 
direct  stiffness  and  finite  difference  techniques  to  simulate  the  3-D  response 


of  granular  array  at  small  strains.  A  more  detailed  discussion  of  this  general 
solution  is  presented  in  Section  4.1. 

The  contact  theory  has  been  repeatedly  used  to  predict  the  elastic  stress- 
strain  properties  of  granular  arrays  of  spheres.  Several  authors  calculated 
the  influence  of  isotropic  confining  pressure  on  Vp  and  Vs  for  various  arrays 
of  smooth  and  rough  spheres,  and  concluded  that  both  velocities  increase  pro¬ 
portionally  to  (aQ)1/6  (Hara,  1935,  Takahashi  and  Sato,  1950,  Gassman,  1951, 
White  and  Sengbush,  1953,  Brandt,  1955).  Of  special  interest  here  are  some 
detailed  analytical  and  experimental  investigations  of  regular  arrays  of  equal 
spheres.  Deresiewicz  (1958)  lists  the  five  stable  regular  arrays  included  in 
Table  1  and  sketched  in  Fig.  11,  which  range  from  the  loosest  simple  cubic, 
(void  ratio,  e  =  0.91)  to  the  densest  pyramidal  (also  called  face  centered 
cubic  array)  and  tethraedral  (also  called  hexagonal  close  packed),  both  with  e 
=  0.35.  More  complete  lists  and  descriptions  of  feasible  regular  arrays  have 
been  presented  by  Filep  (1936),  Brown  (1978)  and  Shahinpoor  (1981).  Table  2 
reproduces  one  of  these  lists  containing  31  arrays,  while  Fig.  12  presents 
elevation  and  plane  views  for  one  of  the  loosest  arrays  of  Table  2  (Cell  No.  2 
with  e  =  1.94).  Deresiewicz  (1958a)  investigated  in  detail  the  simple  cubic 
array  subjected  to  an  initial  isotropic  loading  followed  by  an  arbitrary  stress 
history.  He  found  that  the  array  is  statically  determinate  in  this  case,  pro¬ 
vided  that  the  stress  field  is  uniform,  with  a  one-to-one  correspondence  be¬ 
tween  the  nine  components  of  the  stress  tensor  and  the  nine  independent  com¬ 
ponents  of  the  contact  forces.  Whitman  et  al.  (1964)  studied  a  2-D  version  of 
the  simple  cubic  array  subjected  to  triaxial  and  confined  compression.  Duffy 
and  Mindlin  (1957),  Duffy  (1959)  and  Hendron  (1963)  investigated  the  densest 
arrays  of  Fig.  11  and  Table  1,  which  are  statically  indeterminate,  including 
some  measurements  of  compressional  (rod)  wave  velocity,  VL,  in  stainless  steel 


granular  bars  loaded  Isotropically  (see  Fig.  13).  As  shown  in  the  figure,  the 
measured  values  of  Vl  are  somewhat  smaller  than  predicted,  with  the  difference 
being  greater  at  small  values  of  a0,  and  with  this  difference  increasing  for 
the  low  tolerance  balls;  at  high  pressures  the  measured  approach  the  pre¬ 
dicted  one.  As  a  result,  the  observed  Vl  a(o0)m,  where  m  >  1/6  =  0.167  pre¬ 
dicted  by  the  theory.  This  difference  is  explained  by  Deresiewicz  (1958),  by 
the  small  differences  in  size  between  the  actual  spheres,  which  results  in  the 
number  of  actual,  load-transmitting  contacts  being  smaller  than  predicted; 
thus,  the  array  is  less  stiff  than  calculated.  When  the  tolerance  becomes 
higher  or  the  pressure  increases,  the  number  of  these  actual  contacts  also 
increases  and  approaches  the  theoretical  value,  and  thus  the  measured  velocity 
also  approaches  the  prediction.  As  discussed  in  Section  2,  values  of  m  =  0.25 
>  0.167  have  also  been  measured  for  Vg  in  isotropically  loaded  sands,  most 
probably  due  to  the  same  reason:  an  increase  in  the  number  of  actual  contacts 
as  o0  increases. 

Several  approaches  have  been  used  to  model  the  effect  of  deviations  from 
regularity  in  arrays  of  equal  or  unequal  spheres.  Smith  et  al.  (1929)  pro¬ 
posed  considering  a  random  array  as  formed  by  clusters  of  loose  and  dense  reg¬ 
ular  arrays,  each  present  in  such  proportion  as  to  yield  the  overall  void 
ratio  or  porosity;  this  idea  was  generalized  by  Munro  and  Jowitt  (1974)  and 
Brown  (1978),  who  used  the  concept  of  maximum  entropy  to  find  the  contribution 
of  each  regular  array.  Ko  and  Scott  (1967)  used  a  similar  procedure  to  inves¬ 
tigate  Che  stress-strain  behavior  under  isotropic  compression;  in  this  study, 
"holey"  models  were  used  in  which  some  of  the  spheres  in  both  the  loose  and 
the  dense  regular  component  array,  are  slightly  smaller  than  the  other  spheres 
In  this  way  the  effect  on  the  bulk  modulus  of  the  increased  number  of  contacts 
caused  by  an  increasing  pressure,  was  incorporated  into  the  model.  Perry  and 


-  VV  p.  N  V-".  -■%  -  w  .%  -  * 


■  «  J  VV  /  ^ 


19 


Brown  (1981)  studied  the  influence  of  having  different  size  spheres  on  the 
coiq>liance  of  the  array.  Davis  and  Deresiewicz  (1977)  investigated  the  com¬ 
pressibility  of  a  3-D  random  array  of  smooth  equal  spheres  subjected  to  iso¬ 
tropic  loading.  Serrano  and  Rodrigues-Ortiz  (1973)  suggested  a  method  for 
generating  random  configurations  of  unequal  disks  or  spheres  having  a  pre¬ 
scribed  grain  size  distribution;  their  work  was  continued  by  the  2-D  numerical 
simulations  by  Cundall  and  by  Oner,  previously  discussed. 

4.1  General  Solution  of  the  Contact  Problem 

As  previously  discussed,  Mindlin  and  Deresiewicz  studied  the  problem  of 
the  load-displacement  behavior  of  two  spheres  in  contact.  In  its  most  general 
formulation,  the  problem  is  to  relate  an  arbitrary  force  time  history  P(t)  = 

A  A  A 

Tx(t)i  +  Ty(t)j  +  N(t)k,  transmitted  through  the  contact  (no  twisting  moments 

* 

are  considered),  to  the  corresponding  displacement  time  history  D(t)  =  6x(t)i 
+  6 y ( t ) j  +  a(t)k  of  center  0  of  one  sphere  relative  to  center  O'  of  the  other 
sphere  (see  Figs.  9  and  14).  Tx,  Ty,  <5X,  5y  are  tangential  forces  and  dis¬ 
placements,  while  N  and  a  are  the  normal  force  and  displacement.  Mindlin  and 
Deresiewicz  established  the  general  framework  for  the  solution  of  this  problem, 
and  they  obtained  closed  form  solutions  for  the  following  particular  cases:  i! 
normal  load  only,  N  t  0,  Tx  =  Ty  =  0  (Hertz  problem);  ii)  N  constant  and  Tx 
increasing  or  decreasing  with  Ty  =  0;  and  iii)  a  normal  load  N  followed  by  an 
oscillating  oblique  force  of  constant  dTx/dN  and  Ty  =  0.  However,  the  general 
problem  of  obtaining  D(t)  for  an  arbitrary,  monotonically  or  cyclically  vary- 
ing  P(t)  was  not  solved  until  very  recently  by  Serldi  and  Dobry  (1984).  This 


solution  of  the  contact  problem,  which  has  been  implemented  by  means  of 
currently  available  computer  code  CONTACT,  is  an  elastic-plastic  incremental 
model,  where  all  yield  surfaces  are  cones  of  angle  f  in  force  space  (N,  Tx, 
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Ty),  see  Fig.  14.  The  copies  translate  without  rotation  and  without  changing 

their  shape  or  size  (kinematic  strain-hardening),  and  the  current  positions  of 

-»■ 

the  apexes  of  the  cones  reflect  the  prior  force  history  P(t).  A  modified  form 

+ 

of  the  normality  rule  is  applicable,  and  the  displacement  increment  dD  is 
computed  as: 


dD 


dN(l-vs  )a 


„  dN  ~ 

Es  k  +  f  57  " 


(9) 


where  Es,  vg,  f  are  the  material  properties  of  the  spheres,  a  =  radius  of 
contact  area  between  the  two  spheres,  H,  HQ,  He  are  tangential  elastic  and 
e lasto-plastic  moduli,  and  dTn,  dTt  are  the  outward  normal  and  tangential  com¬ 
ponents  (with  respect  to  the  yield  surface)  of  the  applied  tangential  force 
*  *  ■* 

increment  dT  =  dTn  +  dTt.  As  demonstrated  by  Seridi  and  Dobry  (1984),  this 

elastic-plastic  general  model  reproduces  identically  all  equations  developed 

by  Mindlin  and  Deresiewicz  for  all  their  particular  cases  of  loading,  and 
+  + 
allows  computing  D(t)  for  any  arbitrary  P(t). 

The  availability  of  this  general  solution  to  the  contact  problem  is  crit¬ 
ical  to  the  mathematical  modelling  of  the  small  strain  stress-strain  response 
of  a  granular  soil  by  models  of  spheres.  The  existing  program  CONTACT  relating 
P  and  D  can  now  be  used  as  a  basic  tool,  in  conjunction  with  direct  stiffness 
or  finite  difference  techniques,  to  study  the  behavior  of  either  regular  or 
random  arrays  of  spheres. 
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Section  5 


DIFFERENTIAL  STRESS-STRAIN  RELATIONS  FOR  REGULAR  ARRAYS  OF  SPHERES 

The  general  solution  of  the  problem  of  the  contact  between  two  spheres  can 
be  used  to  derive  incremental  stress-strain  relationships  for  regular  arrays  of 
spheres.  These  stress-strain  relationships  are  discussed  in  this  section,  with 
particular  emphasis  on  the  behavior  under  isotropic  loading  followed  by  very 
small  but  arbitrary  stress  and  strain  increments,  and  for  the  following  regular 
arrays: 

i)  simple  cubic  array  (sc,  see  Figs.  11a  and  17),  discussed  in  Section  5.1 

ii)  body  centered  cubic  array  (bcc,  see  Fig.  15);  discussed  in  Section  5.2 

iii)  face  centered  cubic  or  pyramidal  array  (fee,  see  Figs,  lid  and  16); 
discussed  in  Section  5.3 

iv)  cubical  tetrahedral  (ct,  see  Fig.  lib)  and  tetragonal-sphenoidal  arrays 
(ts,  see  Fig.  11c);  discussed  in  Section  5.4. 

Most  of  these  incremental  stress-strain  relations  were  taken  from  Deresiewicz 
(1958),  Duffy  and  Mindlin  (1957)  and  Moklhouf  and  Stewart  (1967).  However, 
some  are  new;  in  particular,  the  body  centered  cubic  array  is  discussed  here 
for  the  first  time. 

In  addition  to  the  five  regular  arrays  listed  above,  a  sixth  regular  array 
is  the  hexagonal  closed  packed  or  tetrahedral  array  (hep,  see  Fig.  lie,  see 
also  Deresiewicz,  1958).  Table  1  lists  the  most  important  parameters  of  all 
six  arrays.  A  comparison  of  the  behavior  of  the  sc,  bcc  and  fee  cubic  arrays 


is  presented  in  Section  5.5 


5.1 


Simple  Cubic  Array  (sc) 

The  simple  cubic  array  sketched  in  Figs.  11a  and  17  is  the  simplest  of  all 
regular  arrays  of  equal  spheres.  One  sphere  of  radius  R  represents  the  whole 
array,  and  for  a  uniform  stress  field  this  is  a  statically  determinate  system 
with  a  one-to-cne  correspondence  between  the  array's  stresses  and  the  contact 
forces  (Deresiewicz ,  1958).  If  the  normal  stresses  parallel  to  the  axes  of  the 
array  are  o^  (i  =  1,2,3,  see  Fig.  17),  the  normal  contact  forces  are:  N<  = 
4r2 

Oil/*).  If  the  shear  stresses  parallel  to  the  axes  of  the  array  are  0£ j 
(i,j  =  1,2,3  and  i  *  j),  the  corresponding  tangential  contact  forces  and  T^j  = 
4r2 

Of j •  These  relations  occur  due  to  equilibrium  and  are  independent  of  the 
previous  history  of  stresses.  Therefore,  they  are  also  valid  for  any  stress  and 
force  increments  at  any  stage  during  the  loading,  provided  that  no  gross  sliding 
of  the  contact  has  taken  place  dN^  =  4r2  dT^j  =  4R2  do^j.  Figure  17  il¬ 

lustrates  the  case  in  which  an  anisotropic  state  of  stresses  is  applied  first, 
with  all  =  0,  followed  by  small  arbitrary  increments  da^^  and  da^j. 

A  similar  set  of  simple  relations  is  valid  in  this  case  between  the  array 
strains  and  the  displacements  between  spheres;  these  relations  are  obtained  for 
a  uniform  strain  field  based  on  simple  geometric  considerations.  If  cu  =  normal 
relative  displacement  of  centers  of  the  two  adjacent  spheres  separated  by  con¬ 
tact  i,  and  <Sij,  <5^  =  tangential  relative  displacements  between  the  same  two 
centers,  then:  =  aj/(2R);  and  yjj  =  2e^j  =  ( 5 ^ j  +  <5j^)/(2R),  where  = 

normal  strain,  and  =  engineering  shear  strain  of  the  array.  Again,  all 
these  relations  are  independent  of  the  prior  history  of  strains  and  are  valid 
for  incremental  displacements  and  strains. 


/  it  \ 

'  '  Indicial  tensor  notation  is  not  used  here, 
sum  of  several  terras. 


that  is,  oil  does  not  imply  a 
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The  elastic  stiffnesses  corresponding  to  small  stress  and  strain  incre¬ 
ments  applied  to  the  array  subsequent  to  an  isotropic  stress  state,  on  =  022  = 
023  “  o0,  oij  =  0  are: 


3  1/3  -2/3  2  1/3 

where  Sml  =  S2222  =  S3333  =  [~J  1 1  ~v  s  J  [o0GsJ 

(ID 

3  1/3  2(1 —v s ) 1 / 3  2  1/3 

s1212  =  s1313  =  s2323  =  — (2-Vs) —  [o0GsJ 


Notice  that  the  stiffness  matrix  is  diagonal,  and  therefore  the  Poisson's 
ratio  of  the  array,  v  =  de22/dsn  =  0,  for  "triaxial"  loading  corresponding  to 
increasing  an  and  constant  o22  =  033  =  a0.  The  array  has  a  v  =  0  quite 
independently  of  the  values  of  aQ  and  vs  (see  Fig.  20). 

Note  that,  for  a  given  o0,  Eqs.  10-11  describe  a  linearly  elastic 
anisotropic  medium.  The  necessary  and  sufficient  conditions  for  this  medium  to 
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the  coefficients  become: 


3  5/5  ,  2  1/3 

smi  -  ifj  l1-^)'273  l«nGsJ 


3  !/3  2  1/3 

s2222  =  i~J  l1-vsJ  2/3  l °22Gs J 


3  1/3  2  1/3 

s3333  =  [~J  C 1_vs J  2/3  [033gsJ 


3  1/3  *  d-vs)1/3  *  Gs _ 

V  2-vs  i  1/3  !  1/3 

drJ  +  CttttJ 


i  *  J 


where  G„  and  y_  are  the  shear  modulus  and  Poisson's  Ratio  of  the  material  of 
&  s 

the  spheres. 

The  array  locks  under  "triaxial"  conditions  (an  increasing  with  022  = 
°33  =  o0  =  constant),  while  it  fails  in  pure  shear  (aj2  increasing  with  all 
other  da  =  0).  Because  the  stiffness  matrix,  [S],  in  Eq.  10  is  diagonal,  the 
corresponding  diagonal  compliance  matrix  is  [C]  =  [S]~l,  with  each  compliance 
being  the  reciprocal  of  the  corresponding  diagonal  stiffness  terra.  That  is, 

Cnn  =  1/^iiii  and  ^ijij  =  1/Gijij>  and  (de)  =  [C]  {da}. 

In  this  simple  cubic  array,  and  again  for  the  case  of  isotropic  loading, 
aii  =  ajj  =  akk=  ao>  f°r  P_  and  S-Waves  propagating  along  principal  axis  i  of 
the  array  (and,  for  the  S-wave,  polarized  parallel  to  principal  axis  j),  the 
wave  propagation  velocities  Vp  and  Vg  are  proportional  to  (a0)*^,  and  thus, 
the  corresponding  constrained  and  shear  moduli,  D  =  pVp2  =  da<j/de*£  and, 

Gmax  =  PVS  =  dOij/dYij.  are  both  proportional  to  (a0)G*33.  As  discussed  in 
Section  4,  a  similar  dependence  of  modulus  on  (o0)G*33  is  also  predicted  for 
other  regular  arrays,  while  laboratory  measurements  on  regular  arrays  and 
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soils  indicate  that  Gmax  a(o0)®*^. 

For  this  same  case  of  isotropic  loading  and  for  a  cubic  array  of  quartz 
spheres,  and  if  the  array  is  loaded  in  pure  shear,  do^j  =  doji,  then  a  thres¬ 
hold  shear  strain,  yc  =  4.5  x  10-^.  (a0)^^  is  predicted  (see  Appendix  B). 

This  expression  was  obtained  using  the  properties  of  quartz  listed  in  Table  3 
with  a0  in  psi  and  yt  inches/inch,  at  which  gross  sliding  occurs  at  con¬ 
tacts  i  and  j,  and  there  is  a  tendency  for  a  change  in  the  geometric  arrange¬ 
ment  of  the  spheres.  This  predicted  relation  between  yt  and  a0  for  a  sc  array 
of  quartz  spheres  is  plotted  in  Fig.  21  while  Fig.  22  shows  the  detailed  shear 
stress-strain  plots  up  to  the  threshold  (failure).  In  the  range  of  practical 
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interest  for  soils,  500  psf  <  a0  <  2000  psf,  the  expression  gives  yt  =  10  = 

10~2%t  which  agrees  very  well  with  the  measured  yt  in  sands  as  discussed  in 
Section  2.  Expressions  for  the  secant  modulus  reduction,  G/Gmax,  and  for  the 
damping  ratio  of  the  array,  versus  strain  increment  y  =  dy^j  were  also  ob¬ 
tained  for  a  cubic  array  of  quartz  spheres  (Dobry  et  al. ,  1982),  and  were  com¬ 
pared  with  actual  measurements  in  sands,  with  good  agreement.  The  correspond¬ 
ing  comparison  for  G/Gmax  versus  y  is  reproduced  in  Fig.  18  for  an  assumed 
rc  =  1.5  x  10_2%. 

The  more  general  case  of  anisotropic  loading  of  the  cubic  array,  with 
°11  *  J22  *  °33»  very  interesting,  as  this  model  crudely  simulates  the 
laboratory  measurements  of  Vp  and  Vg  on  anisotropically  loaded  sands  discussed 
in  Section  2.  For  a  P-wave  propagating  parallel  to  the  I-axis  of  the  array, 
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the  predicted  expressions  for  D  =  da^/de^  and  Vp  are: 
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D  =  [(3)1/3/2j[Es/(l-vs)  ](0li) 
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Vp  =■  (D/p)l/2 


(14) 


.  *  M  -J • 


where  Es,  vs  =  elastic  constants  of  the  spheres.  Therefore,  both  D  and  Vp  are 
functions  only  of  the  normal  stress  an  in  the  direction  of  propagation,  and 
do  not  depend  on  the  other  two  array  stresses  a j j  and 

For  an  S-wave  propagating  parallel  to  the  i-axis  of  the  anisotropically 
loaded  array  and  with  motions  polarized  parallel  to  the  j-axis  of  the  array, 
the  corresponding  expressions  for  Gmax  =  da^j/dy^j  and  Vs  are: 


[3(1  — v g ) ]1/3(Es)2/3 
(2-Vo)(l+v„) 


1/3  1/3 

aii  aii 
1/3  1/3 

aii  +  aii 


(Wp)1/2 


Gmax  anc^  vs  are  functions  only  of  the  normal  stresses  in  the  direction  of 
propagation  (an)  and  polarization  (ajj),  and  do  not  depend  on  the  third,  out 
of  plane  array  stress  o^.  Furthermore,  as  Eqs.  15  are  symmetric  with  respect 
to  and  ajj,  the  values  of  Vg  and  Graax  do  not  change  if  the  directions  of 
propagation  and  polarization  are  interchanged. 

These  conclusions  for  the  simple  cubic  array,  that  Vp  depends  only  on  a^^, 
and  Vg  depends  only  on  an  and  Ojj,  are  identical  to  the  experimental  findings 
of  Roesler,  Knox  et  al. ,  Kopperman  et  al. ,  and  Yu  and  Richart  on  anisotropi¬ 
cally  loaded  sands,  previously  discussed  in  Section  2.  The  symmetry  of  an  and 
ojj  in  Eq .  15  is  also  present  as  a  symmetry  of  oa  and  in  empirical  Equation 
4,  obtained  by  Knox  et  al.  from  their  sand  measurements. 

Of  course,  Eqs.  14  and  15  cannot  be  used  directly  for  quantitative  pre¬ 
dictions  in  sands,  as  they  give  D  a  and  Gmax  a  for  the  isotropic 

case,  while  D  and  Gmax  a(a0)^^  in  actual  sands.  It  is  interesting  to  modify 

Eq.  15  to  make  it  consistent  with  this  empirical  fact,  by  replacing  a-n^^, 

1  /I 

o-H  '  by  Uii 


1/2  1/2 

'  ,  Ojj  '  ,  and  then  comparing  measurements  and  predictions.  The 


new  equation  is  Gmax  =  j0»5/  (aa®‘  ^+<jj  ^ )  >  where  M  =  constant.  It 

is  useful  to  specialize  this  expression  for  the  biaxial  case,  cj££  =  on  =  ov, 
Ojj  =  o^k  =  033  =  0h»  to  be  able  to  compare  it  with  the  empirical  Eq.  6  ob¬ 
tained  by  Yu  and  Richart  for  sand.  If  K  =  Oii/ojj  =  oil/o33>  the  new  equation 
becomes:  Gmax  =  Nov^*  ^  [2K®*  ^5/ (1+K^1*  ^ )  ] ,  where  N  =  0.5M.  It  is  conven¬ 

ient  to  define  the  normalized  parameter  G  =  ^^/(Nov^^oh^*^),  where  G  =  1 
for  K  =  1.  The  theoretical  expression  G  =  2kO*25/(i+k0.5)  has  been  plotted  in 
Fig.  19.  The  corresponding  empirical  expression  G  =  1-0.18  obtained  from 

Eq.  6,  has  also  been  superimposed  on  Fig.  19  for  typical  values  I^ax  =  3  and  4, 
The  trends  of  the  predicted  and  measured  curves  are  the  same  in  Fig.  19,  with 
the  laboratory  results  showing  a  somewhat  faster  decrease  in  G  as  Kn  increases. 

The  fact  that  the  crude  particulate  model  used  here  is  capable  of  predic¬ 
ting  the  lack  of  influence  of  the  two  stresses  perpendicular  to  propagation  on 
Vp  (Eq.  14),  and  of  the  out-of-plane  stress  on  Vg  and  Gmax  (Eq.  15),  as  well 
as  the  general  trend  of  the  relationship  between  Gmax  and  the  in-plane  stresses 
(Fig.  19),  is  extremely  encouraging.  The  main  advantage  of  the  cubic  array 
used  here  is  its  simplicity,  but  of  course  this  model  is  still  far  from  repre¬ 
senting  real  sand.  One  deficiency  (which  It  shares  with  other  regular  arrays), 
is  that  in  the  general  case  the  array  itself  is  inherently  anisotropic  even 
when  isotropically  loaded  (crystal-type  behavior);  that  is,  Gmax  and  other 
"elastic”  stress-strain  parameters  are  somewhat  different  for  shear  stresses 
corresponding  to  axes  which  are  different  from  the  structural  axes  (axes  of 
symmetry  of  the  array  1,  2,  3  selected  in  Fig.  17).  However,  when  the  material 
of  the  spheres  has  vg  =  0,  the  array  is  isotropic  when  isotropically  loaded. 
Also,  the  array  locks  when  a  "triaxial"  test  is  conducted  on  it  along  any  of 
its  structural  axes,  Instead  of  yielding  and  eventually  failing  in  shear  as  it 
happens  with  actual  granular  materials. 


5.2  Body  Centered  Cubic  Array  (bcc) 


The  body  centered  cubic  array,  sketched  in  Fig.  15,  was  the  next  regular 
array  studied.  It  is  also  represented  by  one  sphere,  and  the  relations  between 
stress  and  contact  forces  can  be  easily  determined  for  one  uniform  stress  field 
of  interest:  isotropic  loading  followed  by  small  stress  increments.  The  coor¬ 
dination  number  (number  of  contacts/sphere)  is  now  eight  instead  of  six  for  the 
simple  cubic  array,  and,  thus,  the  computations  are  somewhat  more  involved.  A 
procedure  analogous  to  that  used  for  analyzing  the  simple  cubic  array  can  be 
followed,  except  that  it  is  now  easier  to  work  directly  with  the  compliance 
matrix,  C^j^,  [C]  =  [S]-^,  instead  of  the  stiffness  matrix  [S]  used  before  in 
Eq.  10. 

For  the  case  of  isotropic  loading,  followed  by  small  stress  increments, 

[Cj  has  the  following  form: 
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Notice  that,  unlike  Eq.  10  for  the  sc  array,  the  compliance  matrix  [C]  in  Eq.  16 
is  not  diagonal.  However,  it  becomes  diagonal  and  it  corresponds  to  an 
isotropic  elastic  medium  with  v=  0,  if  vs  =  0  similarly  to  that  found  for  the  sc 
in  Section  5.1.  For  an  initial  cross-anisotropic  or  biaxial  loading,  033  =  oQ  + 
oa  and  an  =  033  =  ao»  followed  byarbitrary,  small  stress  increments*,  the  form 
of  [C]  is  still  that  of  Eq.  16,  and  the  compliances  in  Eq.  16  are: 
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See  footnote,  Appendix  A. 3. 
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As  it  can  be  seen  in  the  above  equations  17-22,  Vg  and  Vp  are  again  propor¬ 
tional  to  (o0)l/6  for  che  body  centered  cubic  array,  since  the  corresponding 
moduli  are  proportional  to  (a0)^/3j  this  is  a  characteristic  common  to  all 
regular  cubic  arrays  (Duffy  and  Mindlin,  1957,  Duffy,  1959,  Makhlouf  and 
Stewart,  1967).  However,  again  it  is  possible  to  modify  the  exponents  empiri¬ 
cally,  and  (cr0)l/3  can  be  replaced  by  (a0)^/2  when  measurements  and  predictions 
are  compared. 

Threshold  strain  calculations  were  performed  for  this  body-centered  array 
for  the  case  of  triaxial  loading,  starting  from  an  isotropic  pressure  a0,  for 
which  the  array  yields  and  fails  (as  the  array  tends  to  lock  under  pure  shear 
loading).  Again,  the  threshold  strain,  Yt>  obtained  for  the  array  was  a  func¬ 
tion  of  the  confining  pressure,  and  for  an  array  of  quartz  spheres 

(using  properties  for  quartz  in  Table  3),  yt  =  3.44  x  10-3  with  o0  in 

psi  and  in  inches/inch.  This  gives  slightly  lower  values  of  Yt  than  ob¬ 
tained  from  the  simple  cubic  array  of  quartz  spheres  in  Section  5.1,  Yt  = 

4.53  x  10  3  aQ2/3>  The  p^ot;s  for  these  two  expressions  of  Yt  are  compared  in 
Fig.  21,  while  Fig.  23a)  presents  detailed  axial  stress-strain  curves  up  to  the 
threshold  (failure).  For  the  usual  range  of  values  of  confining  pressure  for 
soils,  both  of  them  agree  well  with  Yt  =  10— 2%  experimentally  observed  in 
sands. 

The  body  centered  cubic  array  also  has  some  deficiencies  when  compared  to 
the  behavior  of  actual  sands.  First,  it  remains  isotropic  even  when  loaded 
under  anisotropic  loads,  as  shown  in  Appendix  A.  Second,  in  this  anisotropic 
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loading  case  the  wave  propagation  velocities  are  not  proportional  to  the 
product  of  the  principal  stresses  in  the  directions  of  propagation  and 
polarization  (as  measured  in  sands),  but  rather  they  are  proportional  to  the 
mean  effective  stress,  as  can  be  verified  from  Eqs.  20-22.  Finally,  this 
array  locks  under  pure  shear  loading,  and  in  fact  it  locks  under  a  number  of 
different  shear  loading  paths  depending  upon  their  orientations  and  initial 
stress  state.  However,  the  bcc  array  adds  to  our  understanding  of  the  general 
problem,  as  it  is  a  medium  dense  (e  =  0.47)  array,  located  within  the  range 
between  the  densest  face  centered  cubic  array  (e  =  0.35)  and  the  loosest 
simple  cubic  array  (e  =  0.91). 

5.3  Face  Centered  Cubic  Array  (fee) 

The  Face  Centered  Cubic  Array  sketched  in  Fig.  16  is  one  of  the  two 
densest  arrays,  and  it  has  been  investigated  by  several  researchers  (Duffy 
and  Mindlin  1957,  Ko  and  Scott  1969,  Hendron  1963).  The  differential  stress- 
strain  relationship  for  this  medium  was  derived  by  Duffy  and  Mindlin  (1957). 
The  array  has  12  contacts  per  sphere,  and  unfortunately  it  is  statically  inde¬ 
terminate  for  most  loading  situations;  as  a  consequence,  closed  form  solutions 
are  available  only  for  the  case  of  isotropic  confining  pressure.  In  the  case 
of  transversely  isotropic  loading  a  qualitative  solution  does  exist,  but  the 
compliances  at  the  contacts  have  not  yet  been  evaluated.  Computation  of  these 
compliances  is  a  formidable  task  due  to  the  indeterminacy  of  the  problem  and 
the  variation  of  the  forces  from  contact  to  contact. 

The  incremental  constitutive  law  under  isotropic  confining  pressure  was 
used  by  Duffy  and  Mindlin  (1957)  to  compare  the  theoretical  and  experimental 
rod  wave  velocities  through  a  bar  composed  of  face  centered  cubic  arrays  of 
spheres  (Fig.  13). 
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increments,  do-n>  doij,  the  stiffness  matrix  has  the  form  shown  below: 
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Similar  to  that  discussed  previously  for  Eq.  16  and  the  bcc  array,  the  stiff" 
ness  matrix  [Sj  in  Eq .  23  becomes  isotropic  and  diagonal  only  if  vg  =  0.  In 
that  case,  all  diagonal  terms  are  identical  and  equal  to: 
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and  the  Poisson's  ratio  of  the  array  is  v  =>  0. 


(27) 


If  Che  face  centered  cubic  array  is  consolidated  under  a  transversely  isotropic 


state  of  stress  with  on  =  o0  +  oa,  022  =  033  =  a0,  and  012  =  013  =  023  =  0* 
the  situation  becomes  more  complicated,  since  in  the  case  of  a  non  isotropic 
loading  the  forces  vary  from  contact  to  contact  and  each  compliance  is  differ¬ 
ent.  To  obtain  a  stress-strain  relation  for  an  anisotropic  loading  path,  the 
derivation  must  be  performed  anew,  distinguishing  between  contacts  with 
different  loading  histories. 

The  differential  stress  strain  law  for  this  case  appears  in  great  detail 
in  the  original  paper  by  Duffy  and  Mindlin  (1957);  it  is  identical  in  form  to 
Eq.  23  ,  except  that  now  the  expressions  for  Sjj^  are  not  known.  Thurston 
(1958)  extended  the  results  of  Duffy  and  Mindlin  to  a  set  of  18  equations  and 
18  unknowns. 

The  fee  array  was  subjected  analytically  to  the  conditions  of  a  triaxial 
test  by  Brauns  (1968)  and  Brauns  and  Leussink  (1970),  who  derived  theoretical 
expressions  between  stress  and  strain  at  finite  levels  for  an  array  of  glass 
spheres  (Fig.  23b).  These  expressions  were  later  compared  to  experimental  data 
obtained  in  triaxial  tests  on  regular  fee  packings  of  glass  and  steel  spheres 
(see  Appendix  B3). 

Thurston  and  Deresiewicz  (1959)  derived  expressions  for  the  uniaxial  com¬ 
pression  of  an  fee  array  when  applied  concurrently  with  a  related  isotropic 
pressure.  Again,  the  theoretical  results  were  compared  with  experimental 
results  obtained  through  compression  of  bars  of  steel  bearing  spheres  arranged 
in  fee  array,  with  good  agreement. 

5.4  Cubical-Tetrahedral  (ct)  and  Tetragonal-Sphenoidal  (ts)  Arrays 

The  elastic  constants  relating  stress  to  strain  increments  for  the  Cubical 
Tetrahedral  and  Tetragonal  Sphenoidal  arrays  that  have  been  consolidated 
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isotropically  were  derived  by  Makhlouf  and  Stewart  (1967).  The  procedure  for 
determining  those  constants  is  the  same  as  in  the  other  arrays  described  in 
detail  by  Duffy  and  Mindlin  (1957). 

The  corresponding  constitutive  law  for  the  Cubical  Tetrahedral  array  has 
the  same  form  as  Eq.  23,  but  now 
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As  we  can  see  from  the  above  equations,  the  Cubical-Tetrahedral  array  differs 
from  the  simple  cubic,  the  body  centered  cubic  and  the  force  centered  cubic 
arrays  in  that  it  does  not  exhibit  cubic  anisotropy,  but  rather  transverse 
or  hexagonal  anisotropy,  as  is  also  the  case  of  the  Hexagonal  Close  Packed 
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array  (Duffy,  1959)  and  of  the  Tetragonal-Sphenoidal  Array. 

The  constitutive  law  for  the  Tetragonal-Sphenoidal  array  appears  also  in 
Makhlouf  and  Stewart  (1967).  However,  not  enough  detail  is  provided  in  this 
reference  for  a  full  understanding  of  the  results,  which  are  quite  complex, 
as  the  representative  unit  prism  is  not  symmetric.  Unfortunately,  the 
original  reference  (Makhlouf,  1963)  could  not  be  found  by  the  authors*,  thus 
preventing  a  better  understanding  of  this  array. 

5.5  Comparison  of  Different  Cubic  Arrays 

The  analytical  results  on  the  three  regular  cubic  arrays  discussed  in 
Sections  5. 1-5. 2-5. 3:  simple  cubic,  body  centered  cubic  and  face  centered 
cubic,  were  compared  as  part  of  the  current  research.  This  was  done  to  gain 
further  insight  into  the  behavior  of  granular  media,  and  as  a  necessary  inter¬ 
mediate  step  toward  the  investigation  of  more  elaborated  and  realistic  partic¬ 
ulate  models. 

All  these  arrays  generally  exhibit  cubic  anisotropy  (crystal-type 
behavior)  under  an  isotropic  confining  pressure  oQ.  In  the  three  arrays,  it 
was  found  that  the  necessary  and  sufficient  condition  for  the  array  to  become 
isotropic  under  aQ  is  for  the  Poisson's  ratio  of  the  spheres,  vs,  to  be  equal 


'■  The  differential  constitutive  laws  for  the  cubical-tetrahedral  and  the  tetra¬ 

gonal-sphenoidal  arrays  are  either  not  applicable  to  our  research  or  they  are 
V  erroneous.  As  one  can  see  from  Eqs.  28-31,  contrary  to  general  belief  (Duffy 

•’  1959),  the  cubical  tetrahedral  array  does  not  become  isotropic  under  isotropic 

loading.  This  is  a  serious  deficiency  vis-a-vis  our  research,  as  sands  are 
isotropic  under  isotropic  load.  Consequently,  the  cubical  tetrahedral  array 
will  not  be  used  here.  As  for  the  tetrogonal  sphenoidal  array,  the  results 
are  not  complete,  and  since  the  representative  volume  element  of  this  array  is 
not  symmetric,  completion  of  the  stress  strain  relation  in  Makhlouf  and 
-  Stewart  (1967)  appears  to  be  a  major  task.  The  inexistence  of  the  primary 

i  source  for  this  array,  Makhlouf  (1963),  made  it  impossible  for  the  authors  to 

clarify  the  above  aspects;  therefore  no  results  will  be  used  here  for  the 
cubical-tetrahedral  and  tetragonal-sphenoidal  array. 
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If  vs  =  0,  the  incremental  stiffness  (and  compliance)  matrix  for  the 


jci 


L 


K 


> 


A* 

v 

t 


V 

V 


V 


$ 


1  is 

!  J 


three  arrays  is  diagonal.  Although  the  Poisson's  Ratio  of  quartz  is  us  a  0.15, 
(Table  3),  is  certainly  different  from  zero,  it  is  low  enough  to  make  this 
"vs  =  0  assumption",  needed  for  isotropy,  a  reasonable  one  for  quartz  sands, 
at  least  as  a  first  approximation.  If  vs  =  0  is  assumed,  the  Poisson's  Ratio 
of  the  array  is  also  computed  to  be  v  =  0  for  the  same  three  arrays.  It  must 
be  note  that  for  the  range  0  _<  vs  <  0.5,  values  up  to  v  =  0.13  are  computed 
for  the  same  arrays  (see  Fig.  20).  Therefore,  the  fact  that  a  value  v  =  0 

results  for  the  array  as  soon  as  vg  =  0  is  assumed  does  not  seem  to  be  so  far 

off  either.  It  is  interesting  to  note  that  measurements  of  Vp  and  Vs  by 
Stokoe  and  his  coworkers  on  an  actual  sand  consolidated  isotropically  also 

provided  a  similarly  low  value  of  v  =  0.10  (Knox  et  al.  1982,  Kopperman  et 

al.  1982,  Lee  1985).  In  any  case,  even  with  vs  *  0,  the  cubic  anisotropy  of 
these  expressions  for  v  arrays  is  not  pronounced;  the  error  resulting  from 
computing  the  moduli  between  the  extreme  values  of  vs  is  smaller  than  3.3% 
(Duffy,  1959). 

The  above  three  cubic  arrays  starting  from  an  isotropic  aQ  state,  were 
loaded  statically  in  triaxial  compression  or  pure  shear  up  to  failure,  that  is 
up  to  gross  sliding,  and  computations  were  performed  and  are  displayed  here  for 
their  stress-strain  curves  and  threshold/failure  strains  (Figs.  21-23).  A 
graph  of  obliquity,  o 22^ ao  at  failure  versus  the  intergranular  friction  coef¬ 
ficient  between  spheres,  f,  was  also  computed  and  is  plotted  in  Fig.  24.  The 
curves  obtained  for  the  arrays  in  this  figure  were  also  compared  to  the 
obliquity  obtained  assuming  the  Mohr-Coulomb  Failure  law: 

a22  _  l+sind>  2,,r  .  <K  .  , 

-  =  - - =  tan  (45  +  X)  With  tandj  =  f 

aQ  l-sin<}>  2  Y 
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To  be  able  to  fail  the  simple  cubic  array  in  triaxial  compression,  this  medium 
was  compressed  by  a  force  parallel  to  one  of  the  face  diagonals  of  the  unit 
volume  of  the  array,  that  is  along  a  [110]  direction  (Deresiewicz ,  1959). 

The  same  cubic  arrays  also  give  excellent  results  when  predicting  the  in¬ 
fluence  of  anisotropic  consolidation  on  shear  wave  velocity;  this  is  shown  in 
Fig.  25  by  a  plot  of  normalized  shear  wave  velocity  vs  stress  ratio  K  =  022/00* 
In  this  plot,  Vg(K)  and  Vg ( 1 )  are  the  values  of  the  shear  wave  velocity,  Vg, 
computed  for  the  anisotropic  case  (K)  and  for  the  isotropic  loading  condition 
(K=l),  respectively,  for  direction  of  propagation  and  polarization  parallel  or 
perpedicular  to  022*  The  same  plot  includes  data  measured  by  Stokoe  et  al. 
(1985)  and  Lee  (1985)  on  dry  sand  in  the  large  cubic  testing  facility  at  the 
University  of  Texas,  with  excellent  agreement  between  the  analytical 
predictions  and  experimental  data. 

The  shear  modulus  at  very  small  strains,  Gmax,  computed  for  these  same 
three  cubic  arrays  under  a  given  isotropic  stress,  a0,  is  plotted  in  Fig.  26 
as  a  function  of  the  coordination  number  (=  number  of  contacts/sphere).  As 
expected,  the  higher  the  coordination  number,  the  stiffer  the  array,  with  es¬ 
sentially  a  linear  relation  between  the  two  parameters;  it  is  interesting  that 
for  a  given  aQ  the  straight  lines  in  Fig.  26  extrapolate  down  to  zero,  suggest¬ 
ing  that  Gmax  is  essentially  proportional  to  the  coordination  number  (a  similar 
plot  appears  in  Yanagisawa,  1983).  Therefore,  adding  contacts  to  the  spheres 
has  the  same  effect  on  the  stiffness  of  the  arrays  as  increasing  the  number  of 
springs  in  a  system  of  equal,  parallel,  elastic  springs.  A  derived  plot  is  the 
graph  between  the  shear  wave  velocity,  Vg  =  (Gmax/ p)1 ^  versus  void  ratio  and 
isotropic  pressure  o0  (Fig.  27a)  for  the  same  three  arrays.  This  last  figure 
is  especially  interesting,  as  the  trend  observed  in  actual  sands  is  very 
similar  (compare  analytical  curves  with  measured  data  in  sands  in  Fig.  27b), 


except  that  the  absolute  values  of  Vs  in  the  real  soils  are  much  smaller,  by  a 
factor  of  two  or  three.  For  example,  for  e  =  0.47,  corresponding  to  the  bcc 
array,  and  cr0  =  30  psi  =  4,320  psf,  Vs  =  1,800  fps  is  predicted  by  the  analyti¬ 
cal  model  in  Fig.  27a  ,  while  Vg  =  1,100  fps  has  been  measured  in  rounded 
grained  sands.  Therefore,  Figs.  26  and  27  strongly  suggest  that  the  dependency 
of  Gmax  and  Vs  on  void  ratio  observed  in  real  soils  is  explained  mainly  by  the 
increase  in  the  number  of  contacts  as  the  void  ratio  decreases. 

Even  though  the  above  results  are  encouraging,  the  regular  arrays  are 
still  very  crude  analytical  models  of  actual  granular  soils,  and  results  such 
as  shown  in  Fig.  27a  are  not  easy  to  interpolate  to  intermediate  void  ratios. 

A  significant  improved  model  is  discussed  in  the  following  section. 


Section  6 


A  MODEL  OF  GRANULAR  SOIL  OF  ARBITRARY  VOID  RATIO 

Smith  et  al.  (1929)  found  experimentally  that  a  random  arrangement  of 
equal  spheres,  after  enough  shaking  and  tapping  has  been  applied  to  it,  seems 
to  be  composed  of  regular  arrays,  representing  dense  and  loose  clusters 
distributed  within  the  random  grandular  medium.  The  measurements  showed  that 
all  spheres  had  between  6  and  12  contacts  per  sphere,  which  corresponds  exactly 
to  the  theoretical  range  for  regular  arrays. 

Additional  experimental  work  by  Bernal  and  Mason  (1960),  Bernal  et  al. 
(1964),  Scott  (1960),  Scott  et  al.  (1964),  Davis  (1974)  and  Shahinpoor  and 
Shahrpass  (1982),  Finney  (1983),  Figs.  28  and  29,  has  confirmed  that  2-D  and 
3-D  random  assemblages  of  equal  spheres  tend  to  crystalize.  Consequently,  at 
the  present  time,  it  is  generally  accepted  that  an  assemblage  of  equal  spheres 
can  be  modelled  by  a  combination  of  regular  arrays,  Finney  (1983),  Backman  et 
al.  (1983). 

In  this  section,  a  model  of  granular  soil  is  proposed  which  consists  of 
clusters  of  the  three  cubic  arrays  discussed  in  Section  5.5,  with  the  addi¬ 
tional  assumption  that  the  spheres  have  vs  =  0.  In  this  model,  the  three  cubic 
arrays,  having  different  void  ratios  and  inherent  stiffnesses,  occur  In  pro¬ 
portions  such  as  to  give  the  desired  macroscopic  void  ratio  of  the  "soil".  In 
Section  6.3,  the  relation  between  Gmax,  void  ratio  and  o0  applied  to  the  "soil" 
is  calculated  using  the  self  consistent  method,  and  the  results  are  compared 
with  Gmax  measured  in  actual  sands. 

6.1  The  Se If  Consistent  Method 

One  of  the  most  commonly  used  procedure  for  describing  the  behavior  of 
macroscopically  isotropic  composite  elastic  media  is  the  "Self  Consistent 


Scheme". 


This  "Self  Consistent  Scheme"  was  first  devised  by  Hershey  (1954)  and 
Kroner  (1958)  as  a  means  to  model  the  behavior  of  isotropic  and  anisotropic 
polycrystalline  materials.  Such  materials  are  just  one  phase  media,  but  be¬ 
cause  of  the  random  or  partially  random  orientation  of  the  crystals.  They  are 
heterogeneous,  the  elastic  properties  vary  with  position  within  the  medium  and 
discontinuities  in  properties  exist  across  some  crystal  interfaces. 

In  these  original  applications  of  the  method  to  polycrystalline  aggregates, 
a  single  anisotropic  crystal  was  viewed  as  a  spherical  or  ellipsoidal  inclusion 
within  an  infinite  medium;  this  infinite  medium  had  the  (still  unknown)  iso¬ 
tropic  elastic  properties  of  the  aggregate.  Then  the  medium,  with  the  inclu¬ 
sion  in  it,  was  subjected  to  a  uniform  stress  or  strain  field  applied  at  large 
distances  from  the  inclusion.  Next,  the  orientation  average  of  the  stress  or 
strain  In  the  inclusion  was  assumed  to  be  equal  to  ("consistent  with")  the 
corresponding  applied  value  of  stress  or  strain.  Thus  the  "self-consistent” 
name  of  the  method.  This  formulation  provided  enough  equations  to  solve  for 
the  isotropic  effective  properties  of  the  medium  (Christensen,  1979). 

Improvement  of  this  self  consistent  scheme  and  its  extension  to  multiphase 
media  are  due  to  Hill  (1965)  and  Budiansky  (1965),  who  developed  the  method  to 
be  used  here.  This  Improved  method  represents  an  approximate  analysis  for  the 
prediction  of  the  overall  (macroscopic)  elastic  moduli  of  a  multiphase  medium 
composed  of  a  coherent  mixture  of  several  isotropic,  linearly  elastic  materials 
or  phases.  The  medium  is  assumed  to  consist  of  contiguous,  irregular  zone  con¬ 
tinuing  these  constituent  materials,  and  the  shapes  of  these  zones  are  assumed 
not  to  deviate  much  from  spherical.  The  spatial  distribution  of  the  phases  Is 
assumed  to  be  such  that  the  composite  medium  is  macroscopically  (i.e.  at  a 
scale  much  larger  than  the  dimensions  of  the  zones),  both  homogeneous  and 


isotropic.  Now,  if  an  N-phase  medium  of  total  volume  V  is  defined,  such  that 
the  aggregate  volume  of  all  zones  containing  the  i^h  phase  is  V-j_,  the  volume 
concentration  is  c^  =  V^/V  and  c^  is  also  equal  to  the  probability  that  any 
arbitrary  point  within  the  medium  is  located  within  a  zone  of  the  material. 

It  should  be  noted  that  in  the  limiting  case  of  very  small  concentrations,  cj , 
c2>  •••  CN-1>  the  first  N-l  phases  will  tend  to  appear  as  isolated  inclusions 
in  a  matrix  consisting  of  the  Nth  phase. 

In  order  to  obtain  to  effective  overall  (macroscopic)  shear,  G*,  and  bulk, 
K*,  moduli  of  the  medium,  a  uniform  stress  field  is  applied  at  its  boundaries. 
Then,  the  stress  and  strain  field,  in  each  of  the  phases  is  evaluated  as  ex¬ 
plained  in  the  next  paragraph.  Once  the  fields  are  determined  for  all  mater- 
ials,  the  effective  moduli,  G  and  K  ,  can  be  calculated  by  equating  the  strain 
energies  of  the  macroscopic  medium  and  of  the  phases.  Again,  the  problem 
reduces  to  a  number  of  coupled  equations  for  K  and  G  ,  which  are  in  terms  of 
the  properties  of  the  individual  materials  and  of  their  volume  concentrations 
(Budiansky,  1965).  This  method  has  been  severely  critized  for  taking  enormous 
liberties  with  the  geometrical  arrangement  of  the  phases  (Christensen,  1979). 

To  calculate  the  elastic  field  in  each  material,  the  geometry  of  the  different 
zones  containing  the  phase  is  successively  rearranged  to  view  the  phase  as  a 
single  inclusion.  However,  the  method  is  relatively  simple  and  in  many  in¬ 
stances,  when  used  with  caution,  gives  very  good  results.  Furthermore,  it  has 
been  proven  by  Hill  (1965)  that  this  "Self  Consistent  Method"  yields  results 
for  G  and  K  which  always  lie  between  the  Voigt  and  Reuss  bounds,  that  is,  the 
spatial  average  of  the  moduli  of  the  phases  (Voigt  bound,  springs-in-parallel) 
and  of  the  reciprocal  of  the  moduli,  or  compliances  of  the  phases  (Reuss  bound, 
springs-in-sands) . 


The  evaluation  of  the  stress  and  strain  fields  in  each  of  the  phases  is 


performed  for  the  isotropic  case  by  the  solution  of  the  problem  of  the  fields 
of  an  ellipsoidal  elastic  inclusion  (Eshelby,  1957).  It  was  shown  by  Eshelby 
that  the  fields  inside  an  ellipsoidal  isotropic  elastic  inclusion  embedded  in 
an  isotropic  elastic  medium  is  uniform;  this  is  an  extremely  important 
conclusion  as  it  eliminates  the  need  for  averaging  the  fields  within  the 
inclusion  phase  and  simplifies  enormously  the  formulation.  Later,  it  was  shown 
that  the  stress  and  strain  fields  inside  an  orthotropic  inclusion  embedded  in 
an  orthotropic  medium  are  also  uniform,  as  long  as  the  cross  section  of  the 
inclusion  is  quadratic  (Kinoshita  and  Mura,  1971). 

The  averaging  of  the  shear  moduli  of  all  phases  by  a  strain  energy  balance 
between  ;he  medium  and  the  inclusions  yields: 
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where  G*,  K*  are  the  desired  macroscopic  moduli;  K^,  (i  =  1,2,...N)  are  the 
moduli  of  the  icn  phase,  c ^  =  V^/V  is  the  volume  concentration,  and  y^,  ev 

i 

are  the  values  of  the  average  shear  strain  and  volumetric  strain  respectively, 

inside  the  phase.  The  parameters  t°  and  o°  are  the  shear  stress  and  isotropic 

o 

pressure  applied  at  the  boundary  of  the  medium. 

The  Elshelby  (1957)  solution  gives 
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where  a  ,3  are  components  of  the  Eshelby  S-tensor;  for  the  case  of  spherical 
inclusions  are: 
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where  v  is  the  macroscopic  Poisson's  ratio  of  the  medium: 
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By  "smearing  out",  that  is,  by  replacing  the  matrix  surrounding  each  inclusion 
(phase)  by  the  desired  resultant  macroscopic  medium,  equations  32  and  33 
reduce  to: 


N 


I 

i=l 


N 


l 

i=l 


ci 


i+a*lF 


-  ij 


i 


(41) 


(42) 


which  are  symmetrical  for  the  various  phases.  Therefore,  Budiansky  (1965)  has 
suggested  to  use  equations  41  and  42  for  arbitrary  concentrations  of  the 
constituents  of  the  composite  medium  as  described  previously.  Furthermore, 
Budiansky  (1965)  simplified  equations  36  and  37  to: 
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A  comparison  of  Eq.  43  with  results  obtained  through  statistical  finite 
element  methods,  suggests  that  the  above  equations  do  indeed  model  the 
continuum  described  previously,  including  the  assumption  that  the  stress  and 
strain  fields  of  the  phases  are  approximately  independent  of  location  (Fig.  30, 
see  Petrakis,  1983). 

Finally,  Equations  38-42  have  to  be  solved  simultaneously  to  yield  the 

if  if 

desired  values  of  the  macroscopic  elastic  moduli,  G  ,  K  .  These  resultant 
macroscopic  moduli  are  estimates  of  the  overall  elastic  constants  of  the 
multiphase  medium,  and,  as  mentioned  before,  they  invariably  lie  between  the 
Reuss  and  Voigt  bounds.  Other  solutions  may  provide  narrower  bounds  for  the 
actual  solution  (Hashin  and  Shtrikman,  1963);  however,  the  Self  Consistent 
solution,  in  certain  cases,  also  falls  between  these  narrower  bounds,  thus 
shwoing  its  capability  for  providing  accurate  results  (Hill  1965). 

6 . 2  The  Model 

The  Self  Consistent  Method  is  applied  here  to  evaluate  the  elastic 
constants  of  a  random  assemblage  of  equal,  rough  elastic  spheres  that  has  been 
consolidated  isotropically  and  has  a  prescribed  mean  void  ratio  e.  The  spheres 
are  assigned  the  elastic  properties  of  quartz,  and  the  assemblage  is  assumed 
to  be  composed  of  random  zones,  with  each  zone  consisting  of  a  large  number  of 
spheres  arranged  in  either  of  three  regular  cubic  arrays. 

Recently,  Shahinpoor  (1981)  modelled  a  random  2-D  array  of  equal  steel 
spheres  as  a  combination  of  Voronoi  cells,  derived  an  expression  for  the 
probability  density  function  of  the  void  ratio,  p(e),  and  checked  experimen¬ 
tally  this  analytical  p(e)  by  means  of  an  optical  scanning  technique,  Fig.  28 


(see  also  Shahinpoor  and  Shahrpass,  1982).  The  expression  for  p(e)  is: 


_ _ X  exp(-Xe) _ 
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X  is  obtained  from  the  mean  void  ratio,  e,  of  the  distribution  p(e).  As  men¬ 
tioned  before,  it  is  reasonable  to  model  a  uniform,  rounded-grained  sand  as  a 
random  combination  of  zones  corresponding  to  regular  cubic  arrays,  and  this 
was  the  approach  taken  in  this  work.  The  sand  medium  is  assumed  to  be  com¬ 
posed  of  regular  arrays  in  the  fashion  of  Figs.  28  and  29,  where  each  randomly 
oriented  Voronoi  polyhedron  is  one  of  these  zones,  and  contains  a  regular  ar¬ 
ray  with  many  spheres.  A  cross  section  of  this  3-D  medium  could  be  visualized 
approximately  by  the  actual  photograph  of  the  2-D  medium  in  Fig.  28b;  in  this, 
the  black  spots  are  spheres  and  the  white  are  voids,  and  zones  of  regular 
packings  can  be  clearly  observed.  The  macroscopic  moduli  of  the  whole  medium 
will  be  determined  from  the  properties  of  these  zones  through  the  self 
consistent  scheme. 

As  a  first  step,  the  probability  density  function  of  the  void  ratio,  p(e), 
Eq.  45,  was  transformed  to  the  probability  density  function  of  the  porosity, 
p(n),  with  the  basic  equation  (Benjamin  and  Cornell,  1970): 


p(n)  =  —  p(e) 


as  the  mean  of  the  porosity  distribution,  n,  coincides  with  the  macroscopic 
(measured)  porosity  of  the  "soil",  unlike  the  mean  void  ratio  e,  which  is  not 
identical  to  the  macroscopic  void  ratio  (Petrakis,  1983). 


Then,  Che  probability  density  function  of  the  porosity,  p(n),  was  dis¬ 
cretized  into  three  segments  of  influence,  corresponding  respectively  to  the 
porosities  of  the  three  regular  cubic  arrays  (see  Fig.  31):  sc  (n  =  0.48), 
bcc  (n  =  0.32)  and  fee  (n  =  0.26).  The  values  of  nm^n  and  nmax  used  for  all 
calculations  were  those  of  the  sc  and  fee  arrays. 

For  example,  for  a  prescribed  macroscopic  porosity  n  =  0.35  correspond¬ 
ing  to  a  mean  void  ratio  e  =  0.54,  the  calculation  illustrated  in  Fig.  31 
allowed  determining  the  following  three  volume  concentrations,  c^: 


Array 

£i 

sc 

0.48 

0.1934 

bcc 

0.32 

0.5921 

fee 

0.26 

0.2143 

The  medium  with  these  three  phases  was  then  subjected  to  an  isotropic 

boundary  confining  pressure,  o°,  and  subsequently  subjected  to  small  boundary 

o 

stress  increments  da-jj,  from  which  the  corresponding  elastic,  very  small 
strain  increments  at  the  boundaries,  and  the  mean  cubic  medium  moduli  K*  and 
G*  were  evaluated. 

If  we  now  assume  that  the  phases  are  quadratic  (elliptical  or  circular) 
in  cross  section  we  can  apply  the  Self  Consistent  method.  The  assumption  that 
the  "zones"  are  quadratic  in  cross-section  is  important,  since  if  the 
resulting  stress  and  strain  fields  are  uniform  (see  Section  6.1)  within  each 
"zone",  the  “zone"  can  be  replaced  by  the  representative  cube  of  each  cubic 
array  (Figs.  15-17)  and  the  corresponding  constitutive  laws  are  given  by  Eqs. 
10,  16  and  23.  Since  in  turn  these  relations  depend  upon  the  pressure  acting 
on  each  inclusion  or  phase,  the  value  of  the  stress  field  at  the  boundary  of 
each  of  these  phases  and  Inside  it  can  be  readily  obtained  from  Eshelby's 
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(1957),  and  Budiansky's  (1965)  results: 


— i  _  oKi  , 
°°  =  00  K* 


1 


Ki 

1  +  1J 


■1 


(47) 


where  the  value  of  a*  depends  upon  the  shape  of  the  zone,  (Eshelby,  1957), 


which  here  has  been  assumed  to  be  spherical  for  simplicity.  Note  that  this 
-i 

value  of  aQ  is  independent  of  the  location  of  the  zone,  and  is  thus  the  same 


for  all  zones  containing  the  same  regular  cubic  array  or  phase. 

Equation  47  is  then  replaced  into  Eqs.  10,  16  and  ,  and  the  problem 
finally  reduces  to  the  solution  of  the  following  equations  for  the  three 
phases: 
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where  =  %(o0),  G^  =  G£(a0)  for  the  three  phases  are  obtained  with  Eqs. 

10,  16  and  23  for  i  =  1,  2,  3,  corresponding,  respectively,  to  the  simple 
cubic,  body  centered  cubic  and  face  centered  cubic  regular  arrays. 

6.3  Application  to  Quartz  Sand 

The  proposed  model  was  evaluated  using  as  input  the  elastic  parameters  of 
quartz  for  the  individual  spheres,  which  are  Es  =  11.0  x  10^  psi  and  vs  =  0.15 
(White,  1965,  Ko  &  Scott,  1967,  Lambe  and  Whitman,  1969,  see  Table  3),  and  for 
a  wide  range  of  isotropic  confining  pressures.  The  values  of  the  computed 
shear  modulus  G*  are  plotted  on  solid  lines  in  Figs.  32  and  33  versus  con¬ 
fining  pressure  aQ  =  o°  for  e  =  0.54  and  e  =  0.46,  respectively.  The  values 

o 

of  the  bulk  modulus,  K  ,  were  also  computed,  and  Fig.  34  contains  a  plot  of 

confining  pressure  versus  volumetric  strain  predicted  by  the  model  for  e  = 

0.54,  where  the  volumetric  strain  was  derived  from  this  computed  bulk  modulus, 
—  o  . 
ev  =  ao/K  • 

As  mentioned  before  in  connection  with  Fig.  27,  regular  arrays  are  much 
stiffer,  up  to  3.5  times  stiffer,  than  actual  uniform,  rounded  sands,  and  this 
constitutes  a  serious  defficiency  of  the  model.  However,  Drnevich  and  Richart 
(1970)  have  suceeded  in  increasing  significantly  the  shear  stiffness  of  dry, 
rounded,  uniform  Ottawa  quartz  sand  by  applying  millions  of  cycles  of  a  shear 
strain  slightly  greater  than  the  threshold  strain  in  a  resonant  column  device. 
Figures  32  and  33  also  include  as  data  points  the  resonant  column  experi¬ 


mental  results  obtained  by  Drenvich  and  Richart,  for  macroscopic  void  ratios, 
e  =  0.54  and  e  =  0.46,  respectively.  The  lower  dotted  line  in  each  figure  cor¬ 
responds  to  virgin  or  uncycled  sand  as  predicted  by  the  Hardin  and  Black  (1966) 


correlation,  which  is  essentially  identical  to  the  Hardin  and  Richart  correla¬ 
tion  depicted  in  Fig.  27(b).  In  both  figures,  two  trends  may  be  clearly 
observed:  a)  as  the  number  of  cycles,  N,  increases, the  test  results  steadily 

approach  the  model  values  until,  at  N  =  22  x  10^  cycles,  in  Fig.  32,  the  agree¬ 
ment  becomes  excellent;  and  b)  the  slope  of  the  line  of  shear  modulus  vs  con¬ 
fining  pressure  decreases  from  about  1/2  in  the  uncycled  state  to  about  1/3 
after,  approximately,  1  x  10^  cycles.  The  reason  for  some  of  the  points  showing 
more  scatter  is  probably  because  those  points  were  cycled  more  than  others  with 
the  same  void  ratio,  or  because  their  void  ratios  were  slightly  differed. 

In  the  above  tests,  the  sand  specimens  were  cycled  at  strains  which, 
although  larger  than  the  threshold  strain,  were  small  enough  so  that  no  sig¬ 
nificant  densification  occurs,  and  indeed  the  change  in  measured  (macroscopic) 
void  ratios  between  virgin  and  cycled  specimens  was  very  little  or  negligible; 
thus,  densification  is  certainly  not  the  explanation  in  the  observed  threefold 
increase  in  the  sand  stiffness.  Drnevich  and  Richart  (1970)  speculated  in 

their  paper  that  the  above  behavior  could  be  due  to  wearing  of  the  contacts, 

increase  of  the  contact  areas  or  formation  of  additional  contacts.  The  authors 
think  that  this  third  reason  explains  the  phenomenom  completely,  as  illustrated 
by  the  comparison  with  the  model  in  Figs.  32-33,  and  by  the  relation  between 
stiffness  and  number  of  contacts  in  Fig.  26.  It  is  known  that  in  a  random 
array  of  spheres  there  may  be  less  contacts  than  in  regular  packings  (see  Smith 
et  al. ,  1929)  and,  furthermore,  it  is  possible  to  have  contacts  which  do  not 
transmit  any  load  (dead  contacts).  By  continuous  cycling  such  as  performed  by 
Drnevich  and  Richart  in  their  tests  these  contacts  were  made  load-transmitting 
and  new  ones  were  formed  until  all  or  most  possible  contacts  were  created  and 
the  stiffness  of  the  sand  coincided  with  that  predicted  by  the  model.  It  is 

interesting  that  the  creation  of  contacts  can  also  be  achieved  by  high 
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isotropic  confining  pressures  (Duffy  and  Mindlin,  1958,  Deresiewicz  1958,  see 
also  Fig.  13). 

The  analytical  model  proposed  herein  describes  exactly  this:  since  it 

assumes  that  the  sand  is  a  random  assemblage  of  regular  arrays  of  spherical 

grains,  it  implies  that  the  number  of  contacts  is  the  maximum  possible. 

Furthermore,  as  shown  by  Duffy  and  Mindlin  (1958)  and  Deresiewicz  (1959),  and 

as  illustrated  by  Fig.  13,  when  the  number  of  contacts  increases  the  pressure 

dependency  of  the  moduli  tends  to  change  from  1/2  to  1/3.  All  of  this  is 

interpreted  by  the  model,  which  is  shown  here  to  represent  the  limiting  state, 

in  terms  of  number  of  contacts  and  of  stiffness,  that  a  sand  can  reach. 

o  _ 

Figure  34  shows  the  confining  pressure  aQ  vs.  volumetric  strain,  ev, 
measured  by  Drnevich  and  Richart  on  the  same  Ottawa  sand  specimens  discussed 
above,  for  a  virgin  specimen  and  for  a  specimen  after  1  x  10^  cycles  of  shear 
strain.  Unfortunately,  no  data  is  available  for  sand  specimens  cycled  with 
more  than  1  x  10^  cycles.  In  the  figure,  lines  have  been  passed  which  approxi¬ 
mately  represent  these  experimental  data.  The  line  predicted  by  the  model  is 
also  plotted,  and  again  it  is  clear  that  the  cycling  is  increasing  the  bulk 
modulus  of  the  soil,  thus  making  it  approach  the  curve  predicted  by  the  model, 
as  the  number  of  contacts  increases  toward  the  theoretical,  maximum  value.  It 
would  be  interesting  to  compare  the  difference  between  the  analytical  results 
and  the  experimental  values  for  both  moduli,  K*  and  G*  at  a  given  cycling 
state.  The  experimental  values  for  the  case  of  hydrostatic  compression  (Fig. 
34)  are  closer  to  the  model  predicted  curve  than  Che  measurements  with  shear 
(Figs. 32-33) ,  for  a  comparable  number  of  cycles.  This  is  explained  by  the 
fact  that,  once  the  spheres  have  approached  each  other  due  to  the  cycles  of 

shearing,  an  increase  in  o°  can  complete  the  formation  of  many  contacts,  thus 

o 

further  increasing  the  stiffness  of  the  soil. 


Section  7 


CONCLUSION 

A  particulate  mechanics  (micromechanical)  model  has  been  developed  for 
describing  the  elastic  response  of  assemblages  of  identical  elastic  spheres  of 
arbitrary  macroscopic  porosity,  n,  subjected  to  an  arbitrary  isotropic 
boundary  pressure,  cs0.  The  model  is  based  on  the  Mindlirr-Deresiewicz  theory 
of  bodies  in  contact  and  takes  into  account  the  spatial  variation  of  porosity. 
The  model  assumes  that  the  assemblage  is  composed  of  random  clusters  of 
several  regular  arrays  of  various  porosities  and  it  computes  the  macroscopic 
moduli  by  means  of  the  Self  Consistent  Method. 

The  predictions  of  the  model,  specialized  for  the  case  of  quartz  spheres, 
were  compared  to  measurements  of  shear  modulus,  Gmax  on  uniform  quartz  sands, 
with  good  qualitative  agreement;  however,  the  analytical  "sands"  were  as  much 
as  3.5  times  stiffer  than  the  actual  soils.  This  is  explained  by  the  fact 
that  sands  have  less  effective  contacts  per  grain  than  theoretically  predicted 
for  a  given  porosity  n.  However,  when  the  sand  is  prestrained  by  millions  of 
shearing  cycles  slightly  above  the  threshold  strain,  Gmax  approaches  the 
theoretical  value  without  changing  n,  as  the  theoretical  number  of  contacts  is 
realized.  Thus,  the  model  exhibits  excellent  agreement  with  results  on 
heavily  prestrained  sand,  and  it  also  provides  upper  bounds  for  small  strain 
shear  and  bulk  moduli  for  all  rounded  uniform  sands. 
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Table  2  Feasible  Regular  Arrays  or  "Cells" 

(Shahinpoor,  1981) 

Coordination  No.  =  N  =  No.  of  Contacts  per  Sphere 
[u,m,St.]  gives  No.  of  Contacts  of  Spheres  with  Layer 
Above,  Same  Layer  and  Layer  Below:  N  =  u  +  m  +  2. 


Young's  Modulus 
Poisson's  Ratio 
Coefficient  of 
Intergranular  Friction 


=  11  x  106 
=  0.15 

=  0.5 


Table  3  Properties  of  Quartz  Used  in  this  Report 

[Note  that  Lambe  and  Whitman  (1969)  reported  a  different  value 
vs  =  0.31  for  quartz,  which  in  turn  was  used  for  the  calculations 
in  Dobry  et  al.  (1982).  The  lower  value  vs  =  0.15  used  herein 
is  more  realistic  and  was  obtained  from  White  (1964)  and  Ko  and 
Scott  (1967)] 
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(Kopperman  et  al.,  1982) 


Stress-Strain  Curves  for  Monotonic  Loading 
of  Dry  Granular  Soils. 
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irucker,  et  al. ;  (b)  Lade  and  Duncan;  (c)  Original  Clam- 
DiMaggio  and  Sandler;  (f)  Zienkiewicz,  et  al. 


Initial  state  of  100-disc 
test:  Isotropic  boundary 
stress 


Tangential  Force  - 
Displacement  Relation  for 
N  Constant,  T  Increasing 
(Dobry  and  Grivas,  1978) 


Velocity,  tt  /se 


Contact  Force  Space  and  Conical  Yield  Surfaces, 
Elastic-Plastic  Incremental  Solution  of  Mindlin' 
Problem  (Seridi  and  Dobry,  1984). 


)  Body-centered  cubic  unit  cell 


b)  Body-centered  cubic  structure 


Fig.  15  (Van  Vlack,  1964) 


)  Face-centered  cubic  unit  cell 


b)  Face-centered  cubic  structure 


Fig.  16  (Van  Vlack,  1964) 
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Fig.  17  Regular  Simple  Cubic  Array 
of  Equal  Spheres 


Secant  Shear  Modulus 
Reduction  Versus  Shear 
Strain:  Comparison  Between 
Calculated  G/Gmax  for  Simple 
Cubic  Array  and  Experimental  m 
Range  for  Sand 

(Qobry  et  al.,  1982)  « 


Experimental  ring*  for  undo 
ISotd  and  Idriss,  1970) 


-Slmplo  cubic  array 
Vt  =1.5  x10l% 
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Simple  Cubic  Array 
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Fig.  19  Influence  of  Stress  Ratio, 

K  =  Oi/O 3  on  Gmax: 

Comparison  Between  Prediction 
From  Simple  Cubic  Array  and 
Experimental  Range  for  Sand 
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Range  for  sands 

G=1-0.l8Kn 

Kn=(K-l)/(Km„-t) 

(Yu  &  Riohart.  1984} 
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Fig.  21  Threshold  Strain  for  Simple  Cubic  (sc)  Body  Centered  Cubic  (bcc) 
and  Face  Centered  Cubic  (fee)  Arrays  of  Quartz  Spheres 


AXIAL  STRAIN  ,  X10  J 


a)  Triaxial  Compression  of  Body  Centered  Cubic  Arrays  of 
Quartz  Spheres 


O  1  2  3  4  5  £l1 

“273 

b)  Triaxial  Compression  of  Face  Centered  Cubic  Arrays  of 
Glass  Spheres  (Brauns  1968) 

Fig.  23  Triaxial  Compression  of  Regular  Arrays 
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g.  28  Histogram  and  Random  T 
Equal  Sized  Spherical 
Shahinpoor  and  Shahrpa 
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Fig.  30  Shear  Strain  Experienced  by  Each  Inclusion,  Yj.> 

as  a  Function  of  its  Shrar  Stiffness,  Oj .  Comparison 
Between  Finite  Element  (FF.M)  and  Self  Consi:  tent 
Method  Results 
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determinate  for  initial  isotropic  and  transversely  isotropic  triaxial  loading, 
only  the  equilibrium  conditions  are  sufficient  for  the  solution  of  this  sub¬ 
problem.  However,  this  case  is  much  more  involved  than  the  simple  cubic  array 
and  tedious  calculations  have  to  be  performed. 

Fig.  A2  shows  one  octant  of  a  sphere  at  the  apex  H  as  well  as  the  point 
of  contact  with  the  “central"  sphere  and  the  applied  and  contact  forces.  This 
octant  of  the  sphere  at  H  will  be  treated  as  the  representative  octant.  Now 
equilibrium  equations  have  to  be  written  for  each  spherical  octant  at  every 
apex  and  the  contact  forces  will  have  to  be  solved  for  each  case  separately. 
For  example  for  apex  H,  (Fig.  A2),  the  equilibrium  equations  are:* 


Z  Fx,  -  o  =>  —  dN33  +  —  dN3!  -  —  dN32  +  7  dPn  +7  dPi3  +  7  dPi2  =  0 

1  VT  VT~  VT~  4  4  4 


(A2) 


Z  Fx  =  0  =>  - -  CIN33  +  -  dN32  +  — ~ —  dP22  +  ~  <^12  +  ~  dF32  ~  0 

2  vr  vr  4  44 


(A3) 


Z  FX-  =  0  =>  dN33 - ~  dN3l - -  <^32  +  ~  ^13  +  7  <8*33  +  7  dP2  3  =  0 

3  vr  vr  vr  4  4  4 


(A4) 


the  solution  of  which  is: 


VT” 

dN31  - - —  (dPu  -  dP 33  +  dP12  ~  dP23> 

o 


(A5) 


,  vr , 

dN32  “  -  (dP^  “  2dP22  +  <^33  "  dPI2  +  2dP 13  -  dP2 3)  (A6) 

24 


9 

dN33 


V3 

-  ( dP 1 1  +  dP22  +  dP33  +  2dPj2  +  2dPi3  +  2dP23) 


(A7) 


*  The  primed  symbols  refer  to  the  local  coordinate  system, 


At  this  point  the  state  of  stress  can  be  defined  and  the  constitutive  law 
may  be  determined  for  each  case  (isotropic  or  cross-anisotropic  triaxial 
confinement) . 


A. 2  Isotropic  State  of  Stress: 


then: 


where  C^jj^  are  Che  expressions  (compliances)  in  the  stress-strain  relations, 


for  example, 

in  eqn. 

(A21), 

(A22) .  In 

matrix 

form 

den 

Cllll 

Cl  122 

c  1 1 33 

0 

0 

0 

dan 

de22 

C2211 
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C2233 

0 

0 

0 

da22 

de33 

C3311 

C3322 

C3333 

0 

0 

0 

da33 

d£12 

0 

0 

0 

c1212 

0 

0 

don 

de13 

0 

0 

0 

0 

C1 313 

0 

dai3 

de23 

0 

0 

0 

0 

0 

c2323 

da23 

(A27) 

In  order  for  the  bcc  array  to  be  isotropic  under  isotropic  loading  the  follow¬ 
ing  condition  must  be  satisfied: 

c1212  =  C1313  =  C2323  =  Can  -  C1122  (A28> 


c2211  =  G1 122  =  C1 133  =  c2233  “  c3311  "  c3322  =  C1 1 22 


(A29) 


The  above  conditions  are  satisfied  only  when  vs  =  0;  furthermore,  in  this 
case  of  the  bcc  array,  as  in  the  sc  array,  Cn22  =  0  and  the  compliance  matrix 
is  diagonal. 

In  this  case: 


1/3 


C1 1 1 1  “  c2222  =  c3333  = 


VT  *V5  Gs2  a0 


(A30) 


Since  the  compliance  matrix  is  diagonal,  its  inverse,  the  stiffness  matrix  is 
easily  computed  by  inverting  each  term:  i.e. 


dcr33 

=  ^^2  0o)1/3 

1 

de33 

dai2 

1 

del  2 

dai3 

1 

dei3 

da23 

1 

de23 

_ 

_  _ 

__  _ 

(A31) 

3  2  1/3 

and  clearly  the  shear  modulus  of  the  array,  G,  Is  G  =  —  (4  3  Gs  aQ)  (A32) 
At  this  point  the  Poisson's  Ratio  of  the  bcc  array  may  be  computed  as  follows 


ubcc 


(A3  3) 


ubcc 


a-,s)2/3 

2_VS 

2(l-vs) 

(1~VS)2/3 

2-vs 

1  /3 

(1-Vg) 

(A34) 


For  different  values  of  vg,  the  Poisson's  Ratio  of  the  bcc  array,  vbcc,  may  bi 
computed;  a  plot  of  V]jCC  versus  the  Poisson's  ratio  of  the  spheres,  vs, 
appears  in  Fig.  (21)  together  with  the  variation  of  vgc  and  vfcc  with  vs  for 


easy  comparison. 


A. 3  Transversely  Isotropic  State  of  Stress  (Triaxial  Loading) 


As  in  the  case  of  the  simple  Cubic  Array,  the  application  of  an  aniso¬ 
tropic  stress  increment  will  result  in  a  variation  of  the  contact  forces  and, 
consequently,  of  the  corresponding  contact  stiffnesses.  Therefore,  in  order 
to  obtain  the  stress  strain  relationships,  the  derivation  must  be  performed 
once  more  distinguishing  between  compliances  at  contacts  with  different 
loading  histories.* 

Consider  now  the  case  of  cross-anisotropic  stress  ("triaxial  test") 
imposed  after  the  array  has  been  subjected  to  an  initial  hydrostatic  stress. 
That  is:  (Fig.  A3) 


at  t  =  tQ  an  =  022  =  033  =  °o 


(A35) 


t  "  ti  on  =  033  =  oc 


(A36) 


022  =  o0  +  oa 


(A3  7) 


The  contact  forces  during  this  Anisotropic  Loading  are: 

,  VT  , 

dN_  =  -  “7“  LdPn  ~  dp33] 

J  1  O 


(A38) 


dN32  =  ~  "T4  fdPH  “  2dh  +  dP33l 


(A39) 


dN33  =  t  ~JJ  “  dPzz  +  dP33] 


(A40) 


in  the  case  of  transversely  isotropic  loading  the  above  equations  simplify  to 


The  computation  of  the  compliances  for  the  case  of  anisotropic  loading  for 
both  the  sc  and  the  bcc  arrays  has  been  done  in  order  for  the  results  to  be 
used  only  for  the  special  case  of  wave  propagation.  This  way,  the  load  was 
assumed  to  reverse  direction,  and  for  this  the  elastic  tangential  compliances 
were  used.  In  the  general  case,  the  load  could  either  increase  or  decrease 
monotonically  and  different  compliances  In  each  case  would  apply. 


dNii  =  0 


,  V6“ 

dN32  =  ±  -  dPa 

12 


,  vr 

dN33  =  ±  dP a 

12 


This  is  Che  case  of  the  two  spheres  in  contact  where  the  normal  force  is 

JL 

increasing  from  NQ  to  NQ  +  N  and  the  tangential  force  from  0  to  T  ,  wit 

B  =  —  >  f  (Mindlin  and  Deresiewicz,  1953) 
dN 

6 .vr>  (  w 

dN  dN33 


(usually  0.5  _<  f  <.  0.8  for  sands) 

in  this  case  the  Normal  Compliance,  Cn,  is: 


r  - - — 

°n  2Gsa 


^  3(l-vg) 

where  a3  =  — -  R(No  +  N) 


and  finally 


q  3(1  — V  e)  .  1  Oa  . 

,3,__^R30or1+i^,i 


/ .  ,2  1/3  -1/3  , 

I -(1~Vs) — ,  r  i  +  I  (£s.)i  I 

4VJG2fln  3  CT°  R 


The  Tangential  Compliance,  Ct,  is: 


CC  ■  4 aS  t'6  H146)[l-(l-9)  ^jyfyy]  } 


where  8  =  — ,  L 
0 


T  ,  T*  T* 

-  and  L  =  - 

fNo  fNo 


However,  in  this  case,  since  the  load  decrement  is  small,  then  L* 
therefore 


2-vs 
Ct  *  4G^ 


in  this  case 


Ct  = 


:~v 


c  2(1 


-v 


i  */3  10  _1/3, 


s  °o 


The  constitutive  law  for  this  particular  loading  may  be  developed 
manner  as  in  the  case  of  the  isotropic  loading,  i.e. 

d<522  =  [“  cn  +  7  Ct]  dP22 
o  3 

finally 


2  i  1  /3  i  oa  ”1/3  ?n  2- 

1e22-  — 1;  *■  )  [«^s>2/3+ - 

3vr  avtgs2o0  3  °0  (1_v 


for  vg  =  0.0 


2  l  1/3  io  “1/3 

de2  2  =  —  [“-T  o — J  [!  +  T  It5-]]  d°22 


V3~  4VTgs2 


s  uo 


3  v°o' 


(A48) 


-  L  +  0, 


(A49) 


(A50) 

in  the  same 


(A51) 


^T]  d022 

(A52) 


(A53) 


in  the  case  of  isotropic  loading  and  vg 


0,  the  above  reduces  to 


VT”  4VT  Gs2aQ 


which  is  the  same  obtained  before,  eqn.  (A30).  Analogously 


den 


=*-t 


VT  A VT  gs2< 


I/3  .  0  -1/3 


2/3 


2-ve 


2(l-vs)1/3 


•]  dc 


(A5! 


for  vs  =  0  this  reduces  to  zero. 


The  Shear  Compliance  is  computed  as  follows: 


dlSij  “  [3  Cn  +  3  Ct;^  dPiJ 


(A5( 


deij 


1/3 


1  ,«a, 


-1/3 


,2/3 


2— v, 


|H  =  W  l4a-r'' 


•] 


(as; 


for  vg  =  0  this  reduces  to 


leij  - -  l 


1  1/3  1  °a  "1/3 

,  d°ij 

3 VT  AVTgs2o0  3  °° 


(A5i 


for  o.  =  0  this  reduces  to  eqn.  (A31).  Therefore,  for  the  case  of  vs  *  0 


2  1  1/3  1  sa  _1/3 

Cl  111  ■  C2222  =  ^3333  =  -  [— ZZ — T — J  [1+T(.T:J]  * 


3\T  4VTGs2a0 


3 


Ci  122  ■  C2233  =  C3311  a  —  l— — — J  [ 1  +  ~  17-JJ 

VT  4VT  Cs2a0  3  °° 


*  [a-vs) 


2/3  _ Z-Vt 


2(l~vs)^/3 


(A60) 


,  1  o-  -1/3 


C1212  s  C1313  -  C2323  =  -  — J  [ 1  ^  ^  I 

3 VT  4V3  Gs2a0  3  0 


(A61) 


In  this  case  the  material  will  be  isotropic  under  cross  anisotropic  loading 
only  if  in  the  compliance  matrix 


C1 1 1 1  "  C1 122  =  C1 2 1 2 


Performing  the  calculations  we  see  that  indeed, 


1  1  I/3  1  oa  ”^3  2/3  2-vs 

C1 111-C1122=  1— — p  J  t1+TlT-J]  + 

3V3  ^  Gs2a0  3  °° 


3  v<V 


(l-vs)l/3 


■]  =  c  12 12 


Therefore,  the  bcc  array  is  isotropic  under  cross  anisotropic  loading  ;  this 
is  a  serious  deficiency  of  the  model  and  should  be  attributed  to  the  symmetry 
of  the  array. 

In  the  case  that  vg  ■  0,  the  array  is  still  isotropic  but  this  time  the 
compliance  matrix  is  again  diagonal: 


for  the  conditions  specified  previously. 


*.  ■  JS  fcSi  A  J»  -*»  .'itW.'i  _V 


APPENDIX  B  SIMULATION  OF  TRIAXIAL  AND  PURE  SHEAR  LOADING  IN  CUBIC  ARRAYS 


The  three  regular  cubic  arrays  will  be  subjected  to  a  finite  loading  in 
such  a  way  so  as  to  cause  failure  when  applied  along  a  principal  direction. 

In  the  case  of  the  simple  cubic  array,  this  corresponds  to  a  pure  shear 
loading,  since  for  a  triaxial  loading  the  array  locks.  Conversely,  the  body 
centered  cubic  and  the  face  centered  cubic  array  will  be  subjected  only  to 
a  triaxial  Loading,  since  for  a  pure  shear  loading  the  array  also  lock. 

In  order  to  obtain  a  finite  stress  strain  relation  in  every  case,  the_ 
co’, i lances  need  to  be  integrated  along  the  loading  path.  Once  this  is  done, 
finite  displacements  are  computed  and  then  the  strains  are  obtained  in  the 
same  manner  as  in  the  infinitesimal  constitutive  laws. 

Finally,  the  load  is  increased  until  failure  (gross  sliding)  occurs  in 
the  array;  In  the  statically  determinate  arrays  gross  sliding  at  a  contact 
translates  into  failure  of  the  cubic  array.  In  the  statically  indeterminate 
face  centered  cubic  (fee)  array,  the  medium  does  not  fail  Immediately,  but 
first  the  number  of  contacts  reduces  from  12  to  8  when  the  fee  array  becomes 
statically  determinate,  and  then  sliding  at  one  contact  becomes  failure. 

B1  The  Simple  Cubic  Array  Subjected  to  Pure  Shear  Loading 

Consider  the  Simple  Cubic  Array  shown  in  Fig.  17  and  consider  a  force  T 
acting  in  the  xj  direction  (Fig.  Al).  The  s.c.  array  is  subjected  to  an 
initial  isotropic  force  N0  and  the  value  of  T  increases  monotonically  from 
zero  to  T*,  where  T*  is  the  value  of  T  that  causes  failure  in  the  array,  while 
Nq  remains  constant.  In  this  case  the  tangential  displacement  is  given  by 
Mindlin  (1949): 


HOMMlU 


fN0  [l  “  (.1  f?a_  J  ] 


(A63) 


whe„  a3 , 


is  the  radius  of  contact. 


Now  y  =  ;  substituting  the  expression  for  a3  into  eqn.  (A63)  and  after  trans 

forming  the  forces  into  stresses  we  obtain: 


3 (2-Ve) 


2  1/3 


s/  *  -  r  uo 


(l-v,)J 


(A64) 


The  above  equation  has  been  plotted  in  Fig.  22  for  different  values  of  a0  and 


f;  obviously  failure  occurs  when 
T  =  f00 

and  at  this  point  the  value  of  the  strain  is* 

2  1/3 

3(2-vs)  *  Oq 

Yf  “  Yt  - - TTr-  - ~) 

(l-v-)1/3  12G-2 


(A65) 


(A66) 


Substituting  the  properties  of  quartz  (Table  3)  in  eqn.  (A66)  we  obtain 


Yt  =  x  10~3  (a0)2/3 
with  aQ  in  psi. 


(A66a) 


B2  The  Body  Centered  Cubic  Array  Subiected  to  Triaxial  Loadin^ 


Consider  the  bcc  array  shown  in  Fig.  15,  and  consider  a  force  Pa  acting 
in  the  X2  direction,  Fig.  Al).  The  bcc  array  is  subjected  to  an  initial 
isotropic  stress  a0  and  the  stress  022  =  °a  *s  increased  monotonically  from  0 

*  The  value  of  vg  =  0.15  used  here  was  obtained  from  White  (1964)  and  other 
sources,  and  is  different  from  vs  =  0.31  used  in  a  previous  report  (Dobry 
et  al. ,  1982).  The  value  vs  »  0.15  is  more  representative  of  quartz;  as  a 
result,  the  values  of  the  threshold  strain  Yf  computed  using  Eq.  A66  and  vs 
3  0.15  are  slightly  different  from  those  originally  obtained  by  Dobry  et  al. 
(1982). 


Co  a  value  oa  at  which  sliding  occurs  in  the  array.  In  this  case,  the  loading 
path  is  as  follows 


at  t 


Pll  ’  p22  =  p33  =  po 


at  t  =  t:  PU  =  P33  =  P0 

p22  -  po  +  pa 

The  contact  forces  are  (eqns,  A5,  A6,  A7 ) 


(A67) 

(A68) 

(A69) 


VT 


dN3i  ° 1  T~  ^dPn  ' <u>33^  =  ° 


(A70) 


VT 


dN32  =  ±  ~24”  r^ll  "  2d§2  +  dP33'j  = 


12 


dPa 


(A71) 


vr 


dN33  =  ±  ~  L dPl  1  +  dP22  +  ^33] 


12 


dPa 


(A72) 


and  the  ratio  between  the  increment  of  the  tangential  force  and  the  increment 
of  the  normal  force,  0,  is 


,  dT  dN32  ^ 


(A73) 


In  this  case  0  >  f,  the  coefficient  of  interparticle  friction,  therefore  the 
values  of  the  compliances  are:  (Jiindlin  and  Deresiewicz,  1953) 
a)  Normal  Compliance,  Cn 
l-vs 


Cn  = 


'n  2Gsa 


(A74) 


this  expression  for  Cn  being  valid  no  matter  what  the  loading  history  of  the 
spheres  is.  Now 

1/3 


a  =»  a0  (1  +  0L) 


(A75) 


where 


3  O  -v  „  ) 


at  n  — . 


3(l-vs)  VT 


0 


n  a.  n 


:>'A 


The  vertical  displacement, 


has  two  components; 


a  6  and  a  6 
33  23 


component: 


2  f  2  t 

S22  -  ^-«33  +  ^«23 


the  6  component  is  computed  as  follows: 


dd33  1-Vo  -1/3 

Cn  =-#=  2G^7<1  + 


d|S33  =  1-vs  f  0N_ 
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1  2gT"[1+I"]  dN 
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(A77) 


(A78) 


(A79) 


(A80) 


(A81) 


The  6^  component  will  be  computed  from  the  tangential  compliance  since  it  is 
a  tangential  displacement;  the  tangential  compliance  is  (Mindlin  and 


Deresiewicz,  1953): 


■»  ’ 


Simplifying: 
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where 


'Co  4Gsao 
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f inally 
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Simplifying  further  and  expressing  the  displacement  in  terms  of  stresses  we 


obtain: 

i  cr  2-vs  r  Oq  VT  °a  3//3  f  V2*  °a 
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In  order  Co  compute  the  vertical  displacement  of  the  array,  622,  we  have  to 
substitute  equations  (A81)  and  (A88)  into  eqn.  (A77): 


2  f  2  | 
622  =  —  633  +  —  623 

^  vr 


6  22  =  4R(1-vs)2/3  [-  -U-°-  -J  [t1+^^J 
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From  the  above  equation,  A90,  the  stress-strain  relation  may  be  computed 
for  the  bcc  array  under  isotropic  loading,  and  it  is: 


1/3 


e22  3  L 


°o 


2/3 


2-vc 


4VT  Gg2 


S  3  <V  J  4  (1_v  } 1 / 3 


vt aa,2/3  .  ,f  ,,  vroa  2/3 


[  [ ( 1  +  ~  ~) 
36  °o 


-  1]  [(1  +  (i  -  1)  — -2-)  -  1]]} 

6  3f  0o 


(A91) 


The  above  equation  is  plotted  in  Fig.  23  for  various  values  of  o0  and  f.  At 
this  point,  the  value  of  oa/o0  which  causes  failure  in  the  array  must  be 
determined.  Failure  is  defined  as  sliding  at  the  contacts;  this  time  since 
the  bcc  array  is  statically  determinate  failure  in  one  contact  implies  failure 
of  the  array.  Furthermore,  because  of  the  symmetry  of  the  array,  failure  will 
occur  simultaneously  at  all  contacts.  Sliding  will  occur  when 
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where  A  is  the  area  of  the  face  of  the  bcc  array.  Then? 


T  =  f(N0  +  N) 


VT”  C.VT~  vT“  . 

-A„22-f(-A3a+TA  ,„) 


VT  A  . 

—  o22  =  f(*3  aa  +  0o) 


°22  ,  VT_  f ,  f 

°o  L  3  3 


and  finally,  oa/o0  at  failure,  (oa/o0)^  is 
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If  02^  =  °o  +  °a  (Cotal  stress),  then 
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In  terms  of  total  stress  (Fig.  24) 

In  order  to  compute  the  strain  at  failure,  e22 
equation  for  (a 22/ao^>  e9n*  (A100)  into  the  stress 
(A91).  Doing  this  we  obtain 
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,  we  must  substitute  the 
strain  relation,  eqn. 


e22=[-fr—  J  li  +— j 

Zf  4YT  Gg2  v^f  4  U 

f  2/3  ,  wy-  2/3 

*[[(1+  — — J  -1]-[1  +  [-  -  1)  -~^-J  -1]]} 

Vz-f  e  vT-f 


(A102) 


for  the  properties  of  quartz  (Lamb  and  Whitman,  1964,  Ko  and  Scott,  1967, 
White,  1964) 


Gc  =  4.783  x  106 


vs  =  0.15 


f  =  0.5 


tua  strain  at  failure  is 


e22^  =  3.438  x  10“3  0Q2/ 3  (<n  percent) 
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and  =  e22f  f°r  vs  = 

which  is  the  threshold  strain  yt  for  the  array.  Thus  an  equation  is  plotted 
in  Pig.  21  together  with  the  other  expressions  for  yt  for  the  other  arrays 
for  easy  comparison. 


B3  The  Face  Centered  Cubic  Array  Subiected  to  Triaxial  Loadin 


The  triaxial  loading  of  a  Face  Centered  Cubic  Array  was  solved  by  Brauns 
and  Leussink  (1970).  In  this  work,  the  stress-strain  relationship  is  not  ob¬ 
tained  for  the  whole  range  of  values  of  o22/o0  but  only  for  those  which  make 
the  array  statically  determinate.  It  is  extremely  hard  to  determine  the 
values  of  compliances  for  cross  anisotropic  loading  (Duffy  and  Mindlin,  1957) 
therefore  once  sliding  occurs  and  the  number  of  contacts  decreases  from  12  to 
8,  the  array  becomes  statically  determinate  and  it  is  possible  to  compute  a 
stress-strain  relation  up  to  the  point  that  the  array  fails  (range  b  -  Fig. 


B2).  Once  Che  array  has  failed,  simple  geometric  considerations  make  the 
computation  over  the  strain  range  c  possible  (Fig.  B2). 


Consider  the  fee  array  in  Fig.  Bla  and  the  cross  anisotropic  loading  as 
shown.  The  free  body  diagram  of  the  octant  of  sphere  A  is  shown  in  Fig.  Bib. 
The  equilibrium  conditions  yield 

N  +  T  =  V2~ R2  01  (A104) 

N  -  T  +  Nj  =  iVT  R2  03  (A  1 0  5 ) 

Also  the  sum  of  the  displacements  around  a  closed  path  must  vanish  (Duffy  and 
Mindlin,  1957),  which  yields 

ai  -  a  +  6  =  0  (A106) 


The  normal  compliance  is  found  to  be 
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and  the  Tangential  Compliance 
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The  strain  is 

en  *  ~  (a  +  6)  (A109) 

2R 

Integrating  the  compliances  and  substituting  the  results  into  eqn.  (A109),  we 
find  after  transforming  the  forces  into  stresses  that  the  strain,  ejj,  is 
given  by 
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this  expression  being  valid  only  for  the  range  of  (on/o33)  in  which  the  array 
is  statically  determinate,  that  is  8  contacts  per  sphere  (range  b  in  Fig.  Bl). 

Failure  is  defined  again  as  sliding,  but  this  time  when  the  whole  array 
fails,  that  is  when  the  number  of  contacts  from  8  reduces  to  none. 

Using  the  same  criteria  as  in  the  other  arrays,  the  critical  stress  ratio 
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at  which  the  array  fails  is 
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The  above  equation  is  plotted  in  Fig.  24. 
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Fig.  B1  a)  Representative  Unit  Volume  of  an  fee  Array  with  State 
of  Stress 

b)  Sphere  Centered  at  Apex  A  with  Applied  Stresses,  Contact 
Forces  and  Displacements  (Brauns  and  Leussink,  1970) 
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Fig.  B2  Triaxial  Compression  of  an  fee  Array  of  Glass  Spheres 
Analytical  and  Experimental  Results 
(Brauns  and  Leussink,  1970) 


